In [2]:
from scipy.linalg import solve
# Define the equation parameters
from scipy.special import erfc
import pandas as pd
import numpy as np
import scipy.integrate
from numpy import exp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.animation as FuncAnimation
from IPython.display import HTML
plt.rcParams['figure.figsize'] = [10, 5]
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
import seaborn as sns
from matplotlib.animation import PillowWriter
In [4]:
# Analytical - du/dt = -v*du/dx+D*d2u/dx2 - c*u +d

Ci = 0; 
Co = 1;
# dCinf/dx = 0;
R = 1; 
c = 0.0
D = 0.01
v = 1.0
d = 0.0
nx = 101;
nt = 1001;

L = 1.0;
T = 1.0;

x = np.linspace(0, L, nx);
t = np.linspace(0, T, nt);


C = np.zeros((nt, nx))

for i in range(nt):  
    for j in range(nx):
        C[i, j] = (0.025/(np.sqrt(0.000625 + 0.02*t[i])))*np.exp((-(x[j] + 0.5 - t[i])**2)/(0.00125 + 0.04*t[i]))
        
        # An = np.exp(-c*t[i]/R)*(1-(erfc((R*x[j]-v*t[i])/(2*np.sqrt(D*R*t[i]))))*0.5 -  0.5*(np.exp(v*x[j]/(D)))*(erfc((R*x[j]+v*t[i])/(2*np.sqrt(D*R*t[i])))))       
        # un = v*(np.sqrt(1+ (4*c*D)/(v**2)))
        # Bn = 0.5*np.exp(((v-un)*x[j])/(2*D))*(erfc((R*x[j]-un*t[i])/(2*np.sqrt(D*R*t[i])))) + 0.5*np.exp(((v+un)*x[j])/(2*D))*(erfc((R*x[j] + un*t[i])/(2*np.sqrt(D*R*t[i]))))
        # print(np.shape(An), np.shape(Bn))
        # C[i, j] = d/c + (Ci - d/c)*An + (Co-d/c)*Bn

U_exact = np.transpose(C);

fig = plt.figure()  
plt.imshow(np.transpose(C), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
plt.colorbar()
plt.title("Exact")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDI.png', dpi = 900)
plt.show()
[1.38389653e-87 1.45844014e-84 1.00924437e-81 ... 4.21713645e-04
 4.13846376e-04 4.06120826e-04]



[0.         0.         0.         ... 0.00038267 0.00039424 0.00040612]
In [363]:
from mpl_toolkits.mplot3d import Axes3D
sns.set(style = "darkgrid")
ax = fig.add_subplot(111, projection = '3d')
X, Y = np.meshgrid(x, np.transpose(t))
sns.set_style("whitegrid", {'axes.grid' : False})
fig = plt.figure(figsize=(12,12))
ax = plt.axes(projection='3d')
ax.scatter(X, Y, np.transpose(U_exact), s=4, c=U_exact, cmap = 'viridis', alpha = 1)
ax.set_xlabel('Position')
ax.set_ylabel('Time')
ax.set_zlabel('U_exact')
plt.show()
$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x}$$$$ 0<x<1 \ \ \ 0<t \leq 1 $$$$ IC: U(x, 0) = \exp(- \frac {(x+0.5)^2}{0.00125}) $$

$$ BC - Left: U(0, t) = \frac{0.025}{\sqrt{0.000625 + 0.02t}} \exp\left[ - \frac{(0.5-t)^2}{(0.00125 + 0.04t)} \right] $$ $$ BC - Right: U(1, t) = \frac{0.025}{\sqrt{0.000625 + 0.02t}} \exp \left[- \frac{(1.5-t)^2}{(0.00125 + 0.04t)} \right] $$ $$ 0<t \leq 1 $$

Analytical Solution - Para et al. 2022

For v = 1.0 and D = 0.01 $$ U(x, t) = \frac{0.025}{\sqrt{0.000625 + 0.02t}} \exp \left[- \frac{(x+0.5-t)^2}{(0.00125 + 0.04t)}\right] $$

In [364]:
# Define the grid
a = -1.0;
b = 0.01;
c = 0.0;
d = 0.0;



nx = 101
nt = 1001
L = 1.0
T = 1.0
# Define x & t
x = np.linspace(0, L, nx)
t = np.linspace(0, T, nt)


# Dirichlet Boundary Condition
left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)

fig = plt.figure()  
plt.imshow(np.transpose(C), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
plt.colorbar()
plt.title("Exact")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDI.png', dpi = 900)
plt.show()


# Discretization of space and time
dx = L/(nx-1)
dt = T/(nt-1)
u = np.exp(-((x+0.5)**2)/(0.00125)) #np.zeros(nx);
def numerical_scheme_CDI(u_CDI, dx, dt, l, r):
    alpha = a*dt/dx
    beta = b*dt/dx**2
    A = np.zeros((len(u_CDI), len(u_CDI)))
    bi = d*dt+u_CDI;

    for j in range(1, len(u_CDI)-1):

        A[j, j-1] = -(beta-alpha/2)
        A[j, j] = 1+2*beta-c*dt
        A[j, j+1] = -(beta+alpha/2)

    N = len(u_CDI)-1
    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0
    
    
    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r


    u_CDI = solve(A, bi)

    return u_CDI

def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r):

    frac = 0.1

    le = -np.log(1-l*frac)
    re = -np.log(1-r*frac)
    ke = -np.log(1-u_eCDI*frac)
    N = len(u_eCDI)-1
    #plt.plot(u_eCDI, label='u_eCDI')
    #plt.plot(ke, label='ke')
    A = np.zeros((len(u_eCDI), len(u_eCDI)))
    udterm = -c
    ydterm = d
    vdterm = -a
    wdterm = b
    #plt.legend()
    uk = udterm*dt*np.ones(N+1); # f
    u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
    vk = -vdterm*dt*np.ones(N+1); # a
    v2k = -wdterm*dt*np.ones(N+1); #d
    wk = wdterm*dt*np.ones(N+1); # b


    v2t = np.zeros(N+1)

    v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

    ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
    be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2));
    ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2)

    de = np.zeros(N+1);

    for t in range(0, (N+1)):
        if(t==N):
            de[t] = ke[t];
        else:
            de[t] = ke[t]*(-u2k[t]*np.exp(ke[t])/2+1)+uk[t]+u2k[t]*np.exp(ke[t])

    bi = de;

    for j in range(1, len(u_eCDI)-1):
        A[j, j-1] = ce[j]
        A[j, j] = 1+ae[j]
        A[j, j+1] = be[j]

    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0


    bi[0] = le
    bi[N] = re


    ke = solve(A, bi)
    u_eCDI = ((1-np.exp(-ke))/frac)
    return u_eCDI

def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r):

    # frac = 0.1


    # ke = -np.log(1-u_eCDCN*frac)
    # N = len(u_eCDCN)-1

    # A = np.zeros((len(u_eCDCN), len(u_eCDCN)))
    # udterm = -c
    # ydterm = d
    # vdterm = -a
    # wdterm = b


    # uk = udterm*dt*np.ones(N+1); # f
    # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
    # vk = -vdterm*dt*np.ones(N+1); # a
    # v2k = -wdterm*dt*np.ones(N+1); #d
    # wk = wdterm*dt*np.ones(N+1); # b


    # v2t = np.zeros(N+1)

    # v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

    # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
    # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2));
    # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2))

    # dcn = np.zeros(N+1);
    # for t in range(0, N+1):
    #     if(t==(N)):
    #         dcn[t] = ke[t];
    #     elif(t==0):
    #         dcn[t] = ke[t];
    #     else:
    #         dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


    frac = 0.1

    le = -np.log(1-l*frac)
    re = -np.log(1-r*frac)
    ke = -np.log(1-u_eCDCN*frac)
    N = len(u_eCDCN)-1

    A = np.zeros((len(u_eCDCN), len(u_eCDCN)))

    v2t = np.zeros(N+1)

    v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) 

    bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t)
    acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4)
    ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) -  ((b*dt)/(8*(dx**2)))*(v2t)

    dcn = np.zeros(N+1);
    for t in range(0, N+1):
        if(t==(N)):
            dcn[t] = ke[t];
        elif(t==0):
            dcn[t] = ke[t];
        else:
            dcn[t] = ke[t] - c*dt + ((d*frac + c)*dt*np.exp(ke[t]))*(1-ke[t]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[t+1] - ke[t-1]) + ((b*dt)/(dx**2))*(ke[t+1]+ke[t-1]-2*ke[t]) - ((b*dt)/(4*(dx**2)))*((ke[t+1] - ke[t-1])**2))
            #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


    bi = dcn;

    for j in range(1, N):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = ccn[j]
        A[j, j] = 1+acn[j]
        A[j, j+1] = bcn[j]

    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = le
    bi[N] = re

    # print(A)
    ke = solve(A, bi)
    u_eCDCN = ((1-np.exp(-ke))/frac)
    return u_eCDCN

def numerical_scheme_CDCN(u_CDCN, dx, dt, l, r):
    alpha = a*dt/(4*dx)
    beta = b*dt/(2*dx**2)

    A = np.zeros((len(u_CDCN), len(u_CDCN)))

    bi = np.zeros(len(u_CDCN))
    for j in range(1, len(u_CDCN)-1):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = -(beta - alpha)
        A[j, j] = 1+(2*beta - c*dt/2)
        A[j, j+1] = -(beta + alpha)

        bi[j] = d*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha)

    N = len(u)-1
    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r



    u_CDCN = solve(A, bi)

    return u_CDCN

def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r):
    frac = 0.1
    le = -np.log(1-l*frac)
    re = -np.log(1-r*frac)
    ke = -np.log(1-u_eQI*frac)
    N = len(u_eQI)-1
    #plt.plot(u_eCDI, label='u_eCDI')
    #plt.plot(ke, label='ke')
    A = np.zeros((N+1, N+1))
    udterm = -c
    ydterm = d
    vdterm = -a
    wdterm = b

    uk = udterm*dt*np.ones(N+1); # f
    u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
    vk = -vdterm*dt*np.ones(N+1); # a
    v2k = -wdterm*dt*np.ones(N+1); #d
    wk = wdterm*dt*np.ones(N+1); # b


    # yk = (bia*frac)*dt;#c+e+uk   *np.ones((1, int(N/2)+1)) #d
    v2t = np.zeros(N+1)
    # print(np.shape(v2t))
    v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2])

    v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3])
    # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)])
    # print(np.shape(v2t))
    #v2t = np.concatenate((v2t, v2t[0]))
    #print(np.shape(v2t))


    cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2)
    bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2)
    aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2)
    fq = vk/(8*dx) + v2k*v2t/((8*dx)**2)

    # USING LAST node Q with other side variation
    cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2)
    bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2)
    anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2)
    fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2)

    dq = np.zeros(N+1);
    for t in range(0, N+1):
        if(t==N):
            dq[t] = (ke[t])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2))
        else:
            dq[t] = ke[t]*((-u2k[t]*np.exp(ke[t])/2)+1)+uk[t]+u2k[t]*np.exp(ke[t]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2]


    bi = dq;

    for j in range(1, N-1):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = cq[j]
        A[j, j] = 1+aq[j]
        A[j, j+1] = bq[j]
        A[j, j+2] = fq[j]


    A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one
    A[N-1, N-1] = 1+anq[N-1]
    A[N-1, N] = bnq[N-1]
    A[N-1, N-3] = fnq[N-1]

    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = le
    bi[N] = re


    # print(A)
    ke = solve(A, bi)
    u_eQI = ((1-np.exp(-ke))/frac)
    return u_eQI

def numerical_scheme_QI(u_QI, dx, dt, l, r):
    alpha = a*dt
    beta = b*dt

    A = np.zeros((len(u_QI), len(u_QI)))
    bi = d*dt+u_QI;
    cQ = 3*alpha/(8*dx) - beta/(dx**2)
    bQ = -7*alpha/(8*dx) - beta/(dx**2)
    aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    eQ = 1*alpha/(8*dx)

    cnQ = 7*alpha/(8*dx) - beta/(dx**2)
    bnQ = -3*alpha/(8*dx) - beta/(dx**2)
    anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    enQ = -1*alpha/(8*dx)


    for j in range(1, len(u_QI)-2):

        A[j, j-1] = cQ
        A[j, j] = 1+aQ
        A[j, j+1] = bQ
        A[j, j+2] = eQ

    N = len(u_QI)-1

    A[N-1, N-2] = cnQ  # applied the a>0 terms for this one
    A[N-1, N-1] = 1+anQ
    A[N-1, N] = bnQ
    A[N-1, N-3] = enQ


    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r



    u_QI = solve(A, bi)

    return u_QI


def numerical_scheme_QCN(u_QCN, dx, dt, l, r):
    alpha = a*dt
    beta = b*dt

    A = np.zeros((len(u_QCN), len(u_QCN)))
    bi = d*dt+u_QCN;
    cQ = 3*alpha/(8*dx) - beta/(dx**2)
    bQ = -7*alpha/(8*dx) - beta/(dx**2)
    aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    eQ = 1*alpha/(8*dx)

    cnQ = 7*alpha/(8*dx) - beta/(dx**2)
    bnQ = -3*alpha/(8*dx) - beta/(dx**2)
    anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
    enQ = -1*alpha/(8*dx)


    for j in range(1, len(u_QCN)-2):
        # A[j, j-2] = -d[j]*dt
        A[j, j-1] = cQ/2
        A[j, j] = 1+aQ/2
        A[j, j+1] = bQ/2
        A[j, j+2] = eQ/2
        bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2

    N = len(u_QCN)-1

    A[N-1, N-2] = cnQ/2  # applied the a>0 terms for this one
    A[N-1, N-1] = 1+anQ/2
    A[N-1, N] = bnQ/2
    A[N-1, N-3] = enQ/2
    bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2


    # Boundary condition
    # A[0, 0] = -3
    # A[0, 1] = 4
    # A[0, 2] = -1
    # bi[0] = 0

    # A[N, N] = 3
    # A[N, N-1] = -4
    # A[N, N-2] = 1
    # bi[N] = 0

    A[0,0] = 1.0
    A[N,N] = 1.0

    bi[0] = l
    bi[N] = r



    u_QCN = solve(A, bi)

    return u_QCN


# u_init = u_exact(x, 1.0)

# Apply the numerical scheme

# U_FD = np.zeros((nx, nt))
U_CDI = np.zeros((nx, nt))
U_CDCN = np.zeros((nx, nt))
U_QI = np.zeros((nx, nt))
U_QCN = np.zeros((nx, nt))
U_eCDI = np.zeros((nx, nt))
U_eQI = np.zeros((nx, nt))

U_eCDCN = np.zeros((nx, nt))


# error_2norm_FD_CDI = np.zeros(nt)
error_2norm_CDI_CDCN = np.zeros(nt)
error_2norm_CDCN_QI = np.zeros(nt)
error_2norm_CDCN_eCDCN = np.zeros(nt)
error_2norm_CDCN_QI = np.zeros(nt)
error_2norm_QI_eQI = np.zeros(nt)

error_2norm_CDI_eCDI = np.zeros(nt)
error_2norm_CDI_QI = np.zeros(nt)
error_2norm_QI_QCN = np.zeros(nt)
error_2norm_CDCN_QCN = np.zeros(nt)
error_2norm_CDI_QCN = np.zeros(nt)


error_2norm_eCDCN_eCDI = np.zeros(nt)
error_2norm_eCDCN_eQI = np.zeros(nt)
error_2norm_eCDI_eQI = np.zeros(nt)

error_2norm_Qexact_QI = np.zeros(nt)
error_2norm_Qexact_eQI = np.zeros(nt)

error_2norm_Qexact_CDCN = np.zeros(nt)

error_2norm_Qexact_CDI = np.zeros(nt)
#error_2norm_FD_CDI = np.zeros(nt)


# u_FD = u;
u_CDI=u;
u_CDCN = u;
u_eQI = u;
u_eCDI = u;

u_QI=u;
u_QCN=u;
u_eCDCN=u;

for n in range(nt):
    # print(n)
    # u_FD = numerical_scheme(u_FD, dx, dt)
    # #print(u)
    # #print(str(n*dt)+" "+str(max(u)))
    # U_FD[:, n] = u_FD
    u_CDI = numerical_scheme_CDI(u_CDI, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_CDI[:, n] = u_CDI

    #error = u_FD - u_CDI
    #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)




    u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c, d, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_eCDI[:, n] = u_eCDI

    error = u_CDI - u_eCDI
    error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_CDCN = numerical_scheme_CDCN(u_CDCN, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_CDCN[:, n] = u_CDCN
    error = u_CDI - u_CDCN
    error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c, d, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_eCDCN[:, n] = u_eCDCN
    error = u_CDCN - u_eCDCN
    error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_QI = numerical_scheme_QI(u_QI, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_QI[:, n] = u_QI
    error = u_CDCN - u_QI
    error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_CDI - u_QI
    error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_eQI[:, n] = u_eQI

    error = u_QI - u_eQI
    error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_eCDI - u_eQI
    error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_eCDCN - u_eQI
    error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_eCDCN - u_eCDI
    error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


    u_QCN = numerical_scheme_QCN(u_QCN, dx, dt, left[n], right[n])
    #print(u)
    #print(str(n*dt)+" "+str(max(u)))
    U_QCN[:, n] = u_QCN
    error = u_CDCN - u_QCN
    error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_CDI - u_QCN
    error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    error = u_QI - u_QCN
    error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

    # Comparing with the exact solution
    error = U_exact[:, n] - u_QI
    error_2norm_Qexact_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


    error = U_exact[:, n] - u_CDCN
    error_2norm_Qexact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


    error = U_exact[:, n] - u_CDI
    error_2norm_Qexact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)
    
    error = U_exact[:, n] - u_eQI
    error_2norm_Qexact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

    
plt.rcParams['font.size'] = 18
# fig = plt.figure()  
# plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
# plt.colorbar()
# plt.title("FD")
# plt.xlabel('Time')
# plt.ylabel('Position')

# fig.set_size_inches(30.,18.)
# plt.savefig('U_FD.png', dpi = 900)
# plt.show()

fig = plt.figure()  
plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("CDI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDI.png', dpi = 900)
plt.show()


fig = plt.figure()  
plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("eCDI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_eCDI.png', dpi = 900)
plt.show()


fig = plt.figure() 
plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("CDCN")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDCN.png', dpi = 900)
plt.show()


fig = plt.figure() 
plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("eCDCN")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_CDCN.png', dpi = 900)
plt.show()


fig = plt.figure() 
plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("eQI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_eQI.png', dpi = 900)
plt.show()

fig = plt.figure() 
plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("QI")
plt.xlabel('Time')
plt.ylabel('Position')

#fig.set_size_inches(30.,18.)
# plt.savefig('U_QI.png', dpi = 900)
plt.show()

fig = plt.figure() 
plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
plt.title("QCN")
plt.xlabel('Time')
plt.ylabel('Position')

# fig.set_size_inches(30.,18.)
# plt.savefig('U_QCN.png', dpi = 900)
plt.show()



fig = plt.figure() 
#plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN')
#plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI')
#plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI')
#plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI')
#plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN')

#plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN')
#plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN')
#plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI')
#plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI')
plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI')
plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI')

plt.plot(error_2norm_Qexact_QI, '*', label = 'Exact_QI')


plt.plot(error_2norm_Qexact_CDI, '+', label = 'Exact_CDI')
plt.plot(error_2norm_Qexact_CDCN, 'x', label = 'Exact_CDDCN')

plt.plot(error_2norm_Qexact_eQI, '-', label = 'Exact_eQI')
plt.title("Errors")
plt.legend()
plt.xlabel('Iterations in Time')
plt.ylabel('Error-2-Norm')

#fig.set_size_inches(30.,18.)
# plt.savefig('Error_Norm.png', dpi = 900)
plt.show()

fig = plt.figure() 
# plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI')
#plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN')
#plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI')
#plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI')
#plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN')

#plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI')
#plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN')
#plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN')
#plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI')

#plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI')
#plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI')
#plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI')
plt.plot(np.log(error_2norm_Qexact_QI), '*', label = 'Exact_QI')
plt.plot(np.log(error_2norm_Qexact_CDI), '+', label = 'Exact_CDI')
plt.plot(np.log(error_2norm_Qexact_CDCN), 'x', label = 'Exact_CDCN')
plt.plot(np.log(error_2norm_Qexact_eQI), '-', label = 'Exact_eQI')

plt.legend()
plt.title("Log errors")
plt.xlabel('Iterations in Time')
plt.ylabel('log-Error-2-Norm')

#fig.set_size_inches(30.,18.)
# plt.savefig('Log-ErrorNorm.png', dpi = 900)
plt.show()


#fig = plt.figure() 
#plt.plot(U_exact[:, -1], '-+',label = 'U_exact')
#plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI')
#plt.legend()
#plt.title("Comparison of U_exact and U_CDI")
#plt.xlabel('Length')
#plt.ylabel('U')

#fig.set_size_inches(30.,18.)
# plt.savefig('Log-ErrorNorm.png', dpi = 900)
#plt.show()
In [365]:
plt.plot(U_eQI[:, -100], label='eQI'); plt.plot(U_exact[:, -100], label='Exact'); plt.legend(); plt.title("At 0.1 seconds - Pe = 100")
Out[365]:
Text(0.5, 1.0, 'At 0.1 seconds - Pe = 100')

Adding MMS code with nx -> 20, 40, 80, 160, 240

For equation with IC and BCs:
$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x}$$$$ IC: U(x, 0) = C_i $$$$ BC - Left: U(0, t) = C_o $$$$ BC - Right: \frac {dU(1, t)}{dx} = 0 $$

Notes on MMS:

* Values should be within the domain of the code: TAU or here u/C lies in the range of 0-5 (0-1 for $$ Pe~1 $$) 
In this case, note that whatever manufactured solution is picked must be continuously differentiable through its second spatial derivative and first temporal derivative. In addition, for illustrative purposes, the manufactured solution is chosen so that time will appear in the boundary conditions. Something of this nature can be used as well: 
T(x, t) = Tg0+Tgxcos(ATgxπx)+Tgtcos(ATgtπt)+Tgtxcos(ATgtxπtx)
The chosen equation is: 
$$ U(x,t) = C_o \cos(\pi x) e^{-t}$$

For x and t from 0 to 1: A = 1 and B = pi/2

In [366]:
xx = np.linspace(0, 1, 101)
tt = np.linspace(0, 1, 1001)
uu = np.zeros((101, 1001))
Co = 2;
for i in range(1001):
    for j in range(101):
        uu[j, i] =  Co*np.cos(np.pi*xx[j])*np.exp(-tt[i])
        
        
plt.imshow(uu, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
plt.colorbar()
Out[366]:
<matplotlib.colorbar.Colorbar at 0x210ac871c88>

Implementing that in the equation for the derivative terms are:

$$ \frac {dU}{dt} = -C_o\cos(\pi x) e^{-t}$$$$ \frac {dU}{dx} = -C_o\pi \sin(\pi x) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o\pi^2 \cos(\pi x) e^{-t}$$

The initial and boundary condition terms: Here left is Dirichlet and right is Neumann

$$ U(x, 0) = C_o\cos(\pi x)$$$$ U(0, t) = C_oe^{-t}$$$$ \frac {dU}{dx} (1, t) = -C_o\sin(\pi 1) e^{-t} = 0$$

The new source term will be:

$$ d = \frac {dU}{dt} + v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2}$$$$ d = -C_o\cos(\pi x) e^{-t} - v*C_o\pi \sin(\pi x) e^{-t} + D*C_o\pi^2 \cos(\pi x) e^{-t}$$$$ d = -C_o e^{-t}(\cos(\pi x) + v*\pi \sin(\pi x) - D*\pi^2 \cos(\pi x)) $$

Notes on MMS for transient flows As a first step, to examine temporal order, spatial discretization must be held constant (fixed grid); then, to examine spatial order, temporal discretization must be held constant (fixed time step). Holding either discretization at a constant level introduces a fixed error that applies to all calculations in that group.
The error in such cases is of the form:

$$ DE = DE_temporal + DE_spatial $$

For transient, keeping one of the errors constant, the other can be examined. However, the error we get is sum of constant term, k, with the error, let's say spatial. Hence, to get the the order, the observed order $ \hat{p} $ must be calculated to remove the constant term with three levels instead of the two required for steady flows as:

$$ \hat{p} = \frac {\ln( \frac{||DE_{l+2}|| - ||DE_{l+1}||}{||DE_{l+1}|| - ||DE_l||})}{\ln (r)} $$

nx nt

25 1000
50 1000
100 1000
200 1000
400 1000
100 250
100 500
100 2000
100 4000

$$ d = -C_o e^{-t}(\cos(\pi x) + v*\pi \sin(\pi x) - D*\pi^2 \cos(\pi x)) $$
In [3]:
# Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. 
def numericalcodes(nx, nt): 
    # Define the grid
    a = -1.0;
    b = 0.01;
    c = 0.0;
    Co = 1.0;
    L = 1.0
    T = 1.0
    dm = np.zeros((nt, nx));
    U_exact = np.zeros((nt, nx));
    # Define x & t
    x = np.linspace(0, L, nx)
    t = np.linspace(0, T, nt)
    for i in range(nt):
        for j in range (nx):
            dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - a*np.pi*np.sin(np.pi*x[j]) - b*(np.pi**2)*np.cos(np.pi*x[j]))
    
    # Dirichlet Boundary Condition
    left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    right = np.ones(nt)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    for i in range(nt):
        for j in range (nx):
            U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i])
    fig = plt.figure()  
    plt.imshow(np.transpose(U_exact), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
    plt.colorbar()
    plt.title("Exact")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    # Discretization of space and time
    dx = L/(nx-1)
    dt = T/(nt-1)
    u = Co*np.cos(np.pi*x) #np.zeros(nx);
    def numerical_scheme_CDI(u_CDI, a, b, c, d,  dx, dt, l, r):
        alpha = a*dt/dx
        beta = b*dt/dx**2
        A = np.zeros((len(u_CDI), len(u_CDI)))
        bi = d*dt+u_CDI;

        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta-alpha/2)
            A[j, j] = 1+2*beta-c*dt
            A[j, j+1] = -(beta+alpha/2)

        N = len(u_CDI)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0


        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r


        u_CDI = solve(A, bi)

        return u_CDI

    def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r):

        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDI*frac)
        N = len(u_eCDI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((len(u_eCDI), len(u_eCDI)))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b
        #plt.legend()
        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2));
        ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2)

        de = np.zeros(N+1);

        for xt in range(0, (N+1)):
            if(xt==N):
                de[xt] = ke[xt];
            else:
                de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt])

        bi = de;

        for j in range(1, len(u_eCDI)-1):
            A[j, j-1] = ce[j]
            A[j, j] = 1+ae[j]
            A[j, j+1] = be[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0


        bi[0] = le
        # bi[N] = re


        ke = solve(A, bi)
        u_eCDI = ((1-np.exp(-ke))/frac)
        return u_eCDI

    def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r):

        # frac = 0.1


        # ke = -np.log(1-u_eCDCN*frac)
        # N = len(u_eCDCN)-1

        # A = np.zeros((len(u_eCDCN), len(u_eCDCN)))
        # udterm = -c
        # ydterm = d
        # vdterm = -a
        # wdterm = b


        # uk = udterm*dt*np.ones(N+1); # f
        # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        # vk = -vdterm*dt*np.ones(N+1); # a
        # v2k = -wdterm*dt*np.ones(N+1); #d
        # wk = wdterm*dt*np.ones(N+1); # b


        # v2t = np.zeros(N+1)

        # v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2));
        # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2))

        # dcn = np.zeros(N+1);
        # for t in range(0, N+1):
        #     if(t==(N)):
        #         dcn[t] = ke[t];
        #     elif(t==0):
        #         dcn[t] = ke[t];
        #     else:
        #         dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDCN*frac)
        N = len(u_eCDCN)-1

        A = np.zeros((len(u_eCDCN), len(u_eCDCN)))

        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) 

        bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t)
        acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4)
        ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) -  ((b*dt)/(8*(dx**2)))*(v2t)

        dcn = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==(N)):
                dcn[xt] = ke[xt];
            elif(xt==0):
                dcn[xt] = ke[xt];
            else:
                dcn[xt] = ke[xt] - c*dt + ((d[xt]*frac + c)*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2))
                #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        bi = dcn;

        for j in range(1, N):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = ccn[j]
            A[j, j] = 1+acn[j]
            A[j, j+1] = bcn[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = le
        # bi[N] = re

        # print(A)
        ke = solve(A, bi)
        u_eCDCN = ((1-np.exp(-ke))/frac)
        return u_eCDCN

    def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r):
        alpha = a*dt/(4*dx)
        beta = b*dt/(2*dx**2)

        A = np.zeros((len(u_CDCN), len(u_CDCN)))

        bi = np.zeros(len(u_CDCN))
        for j in range(1, len(u_CDCN)-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = -(beta - alpha)
            A[j, j] = 1+(2*beta - c*dt/2)
            A[j, j+1] = -(beta + alpha)

            bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha)

        N = len(u)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_CDCN = solve(A, bi)

        return u_CDCN

    def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r):
        frac = 0.1
        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eQI*frac)
        N = len(u_eQI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((N+1, N+1))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b

        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        # yk = (bia*frac)*dt;#c+e+uk   *np.ones((1, int(N/2)+1)) #d
        v2t = np.zeros(N+1)
        # print(np.shape(v2t))
        v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2])

        v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3])
        # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)])
        # print(np.shape(v2t))
        #v2t = np.concatenate((v2t, v2t[0]))
        #print(np.shape(v2t))


        cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2)
        bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2)
        aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2)
        fq = vk/(8*dx) + v2k*v2t/((8*dx)**2)

        # USING LAST node Q with other side variation
        cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2)
        bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2)
        anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2)
        fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2)

        dq = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==N):
                dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2))
            else:
                dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2]


        bi = dq;

        for j in range(1, N-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cq[j]
            A[j, j] = 1+aq[j]
            A[j, j+1] = bq[j]
            A[j, j+2] = fq[j]


        A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anq[N-1]
        A[N-1, N] = bnq[N-1]
        A[N-1, N-3] = fnq[N-1]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = le
        # bi[N] = re


        # print(A)
        ke = solve(A, bi)
        u_eQI = ((1-np.exp(-ke))/frac)
        return u_eQI

    def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r):
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QI), len(u_QI)))
        bi = d*dt+u_QI;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QI)-2):

            A[j, j-1] = cQ
            A[j, j] = 1+aQ
            A[j, j+1] = bQ
            A[j, j+2] = eQ

        N = len(u_QI)-1

        A[N-1, N-2] = cnQ  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ
        A[N-1, N] = bnQ
        A[N-1, N-3] = enQ


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_QI = solve(A, bi)

        return u_QI


    def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r):
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QCN), len(u_QCN)))
        bi = d*dt+u_QCN;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QCN)-2):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cQ/2
            A[j, j] = 1+aQ/2
            A[j, j+1] = bQ/2
            A[j, j+2] = eQ/2
            bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2

        N = len(u_QCN)-1

        A[N-1, N-2] = cnQ/2  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ/2
        A[N-1, N] = bnQ/2
        A[N-1, N-3] = enQ/2
        bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_QCN = solve(A, bi)

        return u_QCN


    # u_init = u_exact(x, 1.0)

    # Apply the numerical scheme

    # U_FD = np.zeros((nx, nt))
    U_CDI = np.zeros((nx, nt))
    U_CDCN = np.zeros((nx, nt))
    U_QI = np.zeros((nx, nt))
    U_QCN = np.zeros((nx, nt))
    U_eCDI = np.zeros((nx, nt))
    U_eQI = np.zeros((nx, nt))

    U_eCDCN = np.zeros((nx, nt))


    # error_2norm_FD_CDI = np.zeros(nt)
    error_2norm_CDI_CDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_CDCN_eCDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_QI_eQI = np.zeros(nt)

    error_2norm_CDI_eCDI = np.zeros(nt)
    error_2norm_CDI_QI = np.zeros(nt)
    error_2norm_QI_QCN = np.zeros(nt)
    error_2norm_CDCN_QCN = np.zeros(nt)
    error_2norm_CDI_QCN = np.zeros(nt)


    error_2norm_eCDCN_eCDI = np.zeros(nt)
    error_2norm_eCDCN_eQI = np.zeros(nt)
    error_2norm_eCDI_eQI = np.zeros(nt)

    error_2norm_Exact_QI = np.zeros(nt)
    error_2norm_Exact_eQI = np.zeros(nt)

    error_2norm_Exact_eCDCN = np.zeros(nt)

    error_2norm_Exact_eCDI = np.zeros(nt)
    error_2norm_Exact_CDCN = np.zeros(nt)
    error_2norm_Exact_QCN = np.zeros(nt)
    error_2norm_Exact_CDI = np.zeros(nt)
    #error_2norm_FD_CDI = np.zeros(nt)


    # u_FD = u;
    u_CDI=u;
    u_CDCN = u;
    u_eQI = u;
    u_eCDI = u;

    u_QI=u;
    u_QCN=u;
    u_eCDCN=u;

    for n in range(nt):
        # print(n)
        # u_FD = numerical_scheme(u_FD, dx, dt)
        # #print(u)
        # #print(str(n*dt)+" "+str(max(u)))
        # U_FD[:, n] = u_FD
        u_CDI = numerical_scheme_CDI(u_CDI,a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDI[:, n] = u_CDI

        #error = u_FD - u_CDI
        #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)




        u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDI[:, n] = u_eCDI

        error = u_CDI - u_eCDI
        error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_CDCN = numerical_scheme_CDCN(u_CDCN, a, b, c,  dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDCN[:, n] = u_CDCN
        error = u_CDI - u_CDCN
        error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDCN[:, n] = u_eCDCN
        error = u_CDCN - u_eCDCN
        error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QI = numerical_scheme_QI(u_QI, a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QI[:, n] = u_QI
        error = u_CDCN - u_QI
        error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QI
        error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eQI[:, n] = u_eQI

        error = u_QI - u_eQI
        error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDI - u_eQI
        error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eQI
        error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eCDI
        error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QCN = numerical_scheme_QCN(u_QCN, a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QCN[:, n] = u_QCN
        error = u_CDCN - u_QCN
        error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QCN
        error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_QI - u_QCN
        error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        # Comparing with the exact solution
        error = U_exact[n, :] - u_QI
        error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQI = error_2norm_Exact_QI[n]
        
        error = U_exact[n, :] - u_QCN
        error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQCN = error_2norm_Exact_QCN[n]
        
        error = U_exact[n, :] - u_eCDCN
        error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDCN = error_2norm_Exact_eCDCN[n]

        error = U_exact[n, :] - u_eCDI
        error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDI = error_2norm_Exact_eCDI[n]
        
        error = U_exact[n, :] - u_CDCN
        error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDCN = error_2norm_Exact_CDCN[n]

        error = U_exact[n, :] - u_CDI
        error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDI = error_2norm_Exact_CDI[n]
        
        error = U_exact[n, :] - u_eQI
        error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreQI = error_2norm_Exact_eQI[n]
    # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150)
    errornorms = np.ones(7)
    errornorms[0] = errorCDI
    errornorms[1] = errorCDCN
    errornorms[2] = errorQI
    errornorms[3] = errorQCN
    errornorms[4] = erroreCDI
    errornorms[5] = erroreCDCN
    errornorms[6] = erroreQI
    
    plt.rcParams['font.size'] = 18
    # fig = plt.figure()  
    # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    # plt.colorbar()
    # plt.title("FD")
    # plt.xlabel('Time')
    # plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_FD.png', dpi = 900)
    # plt.show()

    fig = plt.figure()  
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    fig = plt.figure()  
    plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eCDI.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eQI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eQI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_QI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_QCN.png', dpi = 900)
    plt.show()



    fig = plt.figure() 
    #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN')
    #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI')
    #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI')
    #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI')
    #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN')

    #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN')
    #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN')
    #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI')
    #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI')
    # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI')
    # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI')

    plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI')


    plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI')
    plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN')
    plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN')
    
    plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI')
    plt.title("Errors")
    plt.legend()
    plt.xlabel('Iterations in Time')
    plt.ylabel('Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Error_Norm.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI')
    #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN')
    #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI')
    #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI')
    #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN')

    #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI')
    #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN')
    #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN')
    #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI')

    #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI')
    plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI')
    plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI')
    plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN')
    plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN')
    plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI')

    plt.legend()
    plt.title("Log errors")
    plt.xlabel('Iterations in Time')
    plt.ylabel('log-Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    plt.show()


    #fig = plt.figure() 
    #plt.plot(U_exact[:, -1], '-+',label = 'U_exact')
    #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI')
    #plt.legend()
    #plt.title("Comparison of U_exact and U_CDI")
    #plt.xlabel('Length')
    #plt.ylabel('U')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    #plt.show()
    return errornorms
In [4]:
# Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. 
def numericalcodes_dirichlet(nx, nt): 
    # Define the grid
    a = -1.0;
    b = 0.01;
    c = 0.0;
    Co = 2.0;
    L = 1.0
    T = 1.0
    dm = np.zeros((nt, nx));
    U_exact = np.zeros((nt, nx));
    # Define x & t
    x = np.linspace(0, L, nx)
    t = np.linspace(0, T, nt)
    for i in range(nt):
        for j in range (nx):
            dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - a*np.pi*np.sin(np.pi*x[j]) - b*(np.pi**2)*np.cos(np.pi*x[j]))
    
    # Dirichlet Boundary Condition
    left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    right = -Co * np.exp(-t)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    for i in range(nt):
        for j in range (nx):
            U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i])
    fig = plt.figure()  
    plt.imshow(np.transpose(U_exact), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
    plt.colorbar()
    plt.title("Exact")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    # Discretization of space and time
    dx = L/(nx-1)
    dt = T/(nt-1)
    u = Co*np.cos(np.pi*x) #np.zeros(nx);
    def numerical_scheme_CDI(u_CDI, a, b, c, d,  dx, dt, l, r):
        alpha = a*dt/dx
        beta = b*dt/dx**2
        A = np.zeros((len(u_CDI), len(u_CDI)))
        bi = d*dt+u_CDI;

        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta-alpha/2)
            A[j, j] = 1+2*beta-c*dt
            A[j, j+1] = -(beta+alpha/2)

        N = len(u_CDI)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0


        A[0,0] = 1.0
        A[N,N] = 1.0

        bi[0] = l
        bi[N] = r


        u_CDI = solve(A, bi)

        return u_CDI

    def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r):

        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDI*frac)
        N = len(u_eCDI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((len(u_eCDI), len(u_eCDI)))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b
        #plt.legend()
        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2));
        ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2)

        de = np.zeros(N+1);

        for xt in range(0, (N+1)):
            if(xt==N):
                de[xt] = ke[xt];
            else:
                de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt])

        bi = de;

        for j in range(1, len(u_eCDI)-1):
            A[j, j-1] = ce[j]
            A[j, j] = 1+ae[j]
            A[j, j+1] = be[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0

        A[0,0] = 1.0
        A[N,N] = 1.0


        bi[0] = le
        bi[N] = re


        ke = solve(A, bi)
        u_eCDI = ((1-np.exp(-ke))/frac)
        return u_eCDI

    def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r):

        # frac = 0.1


        # ke = -np.log(1-u_eCDCN*frac)
        # N = len(u_eCDCN)-1

        # A = np.zeros((len(u_eCDCN), len(u_eCDCN)))
        # udterm = -c
        # ydterm = d
        # vdterm = -a
        # wdterm = b


        # uk = udterm*dt*np.ones(N+1); # f
        # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        # vk = -vdterm*dt*np.ones(N+1); # a
        # v2k = -wdterm*dt*np.ones(N+1); #d
        # wk = wdterm*dt*np.ones(N+1); # b


        # v2t = np.zeros(N+1)

        # v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2));
        # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2))

        # dcn = np.zeros(N+1);
        # for t in range(0, N+1):
        #     if(t==(N)):
        #         dcn[t] = ke[t];
        #     elif(t==0):
        #         dcn[t] = ke[t];
        #     else:
        #         dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDCN*frac)
        N = len(u_eCDCN)-1

        A = np.zeros((len(u_eCDCN), len(u_eCDCN)))

        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) 

        bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t)
        acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4)
        ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) -  ((b*dt)/(8*(dx**2)))*(v2t)

        dcn = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==(N)):
                dcn[xt] = ke[xt];
            elif(xt==0):
                dcn[xt] = ke[xt];
            else:
                dcn[xt] = ke[xt] - c*dt + ((d[xt]*frac + c)*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2))
                #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        bi = dcn;

        for j in range(1, N):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = ccn[j]
            A[j, j] = 1+acn[j]
            A[j, j+1] = bcn[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0

        A[0,0] = 1.0
        A[N,N] = 1.0

        bi[0] = le
        bi[N] = re

        # print(A)
        ke = solve(A, bi)
        u_eCDCN = ((1-np.exp(-ke))/frac)
        return u_eCDCN

    def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r):
        alpha = a*dt/(4*dx)
        beta = b*dt/(2*dx**2)

        A = np.zeros((len(u_CDCN), len(u_CDCN)))

        bi = np.zeros(len(u_CDCN))
        for j in range(1, len(u_CDCN)-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = -(beta - alpha)
            A[j, j] = 1+(2*beta - c*dt/2)
            A[j, j+1] = -(beta + alpha)

            bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha)

        N = len(u)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0

        A[0,0] = 1.0
        A[N,N] = 1.0

        bi[0] = l
        bi[N] = r
        
        u_CDCN = solve(A, bi)

        return u_CDCN

    def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r):
        frac = 0.1
        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eQI*frac)
        N = len(u_eQI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((N+1, N+1))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b

        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        # yk = (bia*frac)*dt;#c+e+uk   *np.ones((1, int(N/2)+1)) #d
        v2t = np.zeros(N+1)
        # print(np.shape(v2t))
        v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2])

        v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3])
        # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)])
        # print(np.shape(v2t))
        #v2t = np.concatenate((v2t, v2t[0]))
        #print(np.shape(v2t))


        cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2)
        bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2)
        aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2)
        fq = vk/(8*dx) + v2k*v2t/((8*dx)**2)

        # USING LAST node Q with other side variation
        cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2)
        bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2)
        anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2)
        fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2)

        dq = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==N):
                dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2))
            else:
                dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2]


        bi = dq;

        for j in range(1, N-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cq[j]
            A[j, j] = 1+aq[j]
            A[j, j+1] = bq[j]
            A[j, j+2] = fq[j]


        A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anq[N-1]
        A[N-1, N] = bnq[N-1]
        A[N-1, N-3] = fnq[N-1]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0

        A[0,0] = 1.0
        A[N,N] = 1.0

        bi[0] = le
        bi[N] = re


        # print(A)
        ke = solve(A, bi)
        u_eQI = ((1-np.exp(-ke))/frac)
        return u_eQI

    def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r):
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QI), len(u_QI)))
        bi = d*dt+u_QI;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QI)-2):

            A[j, j-1] = cQ
            A[j, j] = 1+aQ
            A[j, j+1] = bQ
            A[j, j+2] = eQ

        N = len(u_QI)-1

        A[N-1, N-2] = cnQ  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ
        A[N-1, N] = bnQ
        A[N-1, N-3] = enQ


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0

        A[0,0] = 1.0
        A[N,N] = 1.0

        bi[0] = l
        bi[N] = r

        u_QI = solve(A, bi)

        return u_QI


    def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r):
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QCN), len(u_QCN)))
        bi = d*dt+u_QCN;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QCN)-2):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cQ/2
            A[j, j] = 1+aQ/2
            A[j, j+1] = bQ/2
            A[j, j+2] = eQ/2
            bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2

        N = len(u_QCN)-1

        A[N-1, N-2] = cnQ/2  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ/2
        A[N-1, N] = bnQ/2
        A[N-1, N-3] = enQ/2
        bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0

        A[0,0] = 1.0
        A[N,N] = 1.0

        bi[0] = l
        bi[N] = r

        u_QCN = solve(A, bi)

        return u_QCN


    # u_init = u_exact(x, 1.0)

    # Apply the numerical scheme

    # U_FD = np.zeros((nx, nt))
    U_CDI = np.zeros((nx, nt))
    U_CDCN = np.zeros((nx, nt))
    U_QI = np.zeros((nx, nt))
    U_QCN = np.zeros((nx, nt))
    U_eCDI = np.zeros((nx, nt))
    U_eQI = np.zeros((nx, nt))

    U_eCDCN = np.zeros((nx, nt))


    # error_2norm_FD_CDI = np.zeros(nt)
    error_2norm_CDI_CDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_CDCN_eCDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_QI_eQI = np.zeros(nt)

    error_2norm_CDI_eCDI = np.zeros(nt)
    error_2norm_CDI_QI = np.zeros(nt)
    error_2norm_QI_QCN = np.zeros(nt)
    error_2norm_CDCN_QCN = np.zeros(nt)
    error_2norm_CDI_QCN = np.zeros(nt)


    error_2norm_eCDCN_eCDI = np.zeros(nt)
    error_2norm_eCDCN_eQI = np.zeros(nt)
    error_2norm_eCDI_eQI = np.zeros(nt)

    error_2norm_Exact_QI = np.zeros(nt)
    error_2norm_Exact_eQI = np.zeros(nt)

    error_2norm_Exact_eCDCN = np.zeros(nt)

    error_2norm_Exact_eCDI = np.zeros(nt)
    error_2norm_Exact_CDCN = np.zeros(nt)
    error_2norm_Exact_QCN = np.zeros(nt)
    error_2norm_Exact_CDI = np.zeros(nt)
    #error_2norm_FD_CDI = np.zeros(nt)


    # u_FD = u;
    u_CDI=u;
    u_CDCN = u;
    u_eQI = u;
    u_eCDI = u;

    u_QI=u;
    u_QCN=u;
    u_eCDCN=u;

    for n in range(nt):
        # print(n)
        # u_FD = numerical_scheme(u_FD, dx, dt)
        # #print(u)
        # #print(str(n*dt)+" "+str(max(u)))
        # U_FD[:, n] = u_FD
        u_CDI = numerical_scheme_CDI(u_CDI,a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDI[:, n] = u_CDI

        #error = u_FD - u_CDI
        #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)




        u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDI[:, n] = u_eCDI

        error = u_CDI - u_eCDI
        error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_CDCN = numerical_scheme_CDCN(u_CDCN, a, b, c,  dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDCN[:, n] = u_CDCN
        error = u_CDI - u_CDCN
        error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDCN[:, n] = u_eCDCN
        error = u_CDCN - u_eCDCN
        error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QI = numerical_scheme_QI(u_QI, a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QI[:, n] = u_QI
        error = u_CDCN - u_QI
        error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QI
        error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eQI[:, n] = u_eQI

        error = u_QI - u_eQI
        error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDI - u_eQI
        error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eQI
        error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eCDI
        error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QCN = numerical_scheme_QCN(u_QCN, a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QCN[:, n] = u_QCN
        error = u_CDCN - u_QCN
        error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QCN
        error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_QI - u_QCN
        error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        # Comparing with the exact solution
        error = U_exact[n, :] - u_QI
        error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQI = error_2norm_Exact_QI[n]
        
        error = U_exact[n, :] - u_QCN
        error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQCN = error_2norm_Exact_QCN[n]
        
        error = U_exact[n, :] - u_eCDCN
        error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDCN = error_2norm_Exact_eCDCN[n]

        error = U_exact[n, :] - u_eCDI
        error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDI = error_2norm_Exact_eCDI[n]
        
        error = U_exact[n, :] - u_CDCN
        error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDCN = error_2norm_Exact_CDCN[n]

        error = U_exact[n, :] - u_CDI
        error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDI = error_2norm_Exact_CDI[n]
        
        error = U_exact[n, :] - u_eQI
        error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreQI = error_2norm_Exact_eQI[n]
    # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150)
    errornorms = np.ones(7)
    errornorms[0] = errorCDI
    errornorms[1] = errorCDCN
    errornorms[2] = errorQI
    errornorms[3] = errorQCN
    errornorms[4] = erroreCDI
    errornorms[5] = erroreCDCN
    errornorms[6] = erroreQI
    
    plt.rcParams['font.size'] = 18
    # fig = plt.figure()  
    # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    # plt.colorbar()
    # plt.title("FD")
    # plt.xlabel('Time')
    # plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_FD.png', dpi = 900)
    # plt.show()

    fig = plt.figure()  
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    fig = plt.figure()  
    plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eCDI.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eQI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eQI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_QI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_QCN.png', dpi = 900)
    plt.show()



    fig = plt.figure() 
    #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN')
    #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI')
    #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI')
    #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI')
    #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN')

    #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN')
    #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN')
    #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI')
    #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI')
    # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI')
    # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI')

    plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI')


    plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI')
    plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN')
    plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN')
    
    plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI')
    plt.title("Errors")
    plt.legend()
    plt.xlabel('Iterations in Time')
    plt.ylabel('Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Error_Norm.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI')
    #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN')
    #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI')
    #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI')
    #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN')

    #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI')
    #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN')
    #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN')
    #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI')

    #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI')
    plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI')
    plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI')
    plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN')
    plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN')
    plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI')

    plt.legend()
    plt.title("Log errors")
    plt.xlabel('Iterations in Time')
    plt.ylabel('log-Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    plt.show()


    #fig = plt.figure() 
    #plt.plot(U_exact[:, -1], '-+',label = 'U_exact')
    #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI')
    #plt.legend()
    #plt.title("Comparison of U_exact and U_CDI")
    #plt.xlabel('Length')
    #plt.ylabel('U')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    #plt.show()
    return errornorms
In [8]:
# Spatial MMS tests
errorsx = np.ones((9, 7))
errornorms = np.ones(7)
nx = 25;
nt = 100;
nxs = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400]
nts = [250, 500, 1000, 2000, 4000]
for i in range(9):
    errornorms = numericalcodes(nx, nt)
    errorsx[i, 0] = errornorms[0]
    errorsx[i, 1] = errornorms[1]
    errorsx[i, 2] = errornorms[2]
    errorsx[i, 3] = errornorms[3]
    errorsx[i, 4] = errornorms[4]
    errorsx[i, 5] = errornorms[5]
    errorsx[i, 6] = errornorms[6]
    nx = nx*2;
print(errorsx)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI')
plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN')
plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI')
plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN')
plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI')
plt.plot(nxs, np.log(errorsx[:, 5]), '-v', label = 'eCDCN')
plt.plot(nxs, np.log(errorsx[:, 6]), '--', label = 'eQI')
plt.xlabel('Grid Numbers')
plt.ylabel('Log Error (MSE)')
plt.legend()

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
px = np.zeros((7, 7))
for g in range(7):
    for k in range(7):
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))
plt.plot(nxs[2:9], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[2:9], px[:, 1], '-.', label = 'CDCN')
plt.plot(nxs[2:9], px[:, 2], '-+', label = 'QI')
plt.plot(nxs[2:9], px[:, 3], '-*', label = 'QCN')
plt.plot(nxs[2:9], px[:, 4], '-x', label = 'eCDI')
plt.plot(nxs[2:9], px[:, 5], '-v', label = 'eCDCN')
plt.plot(nxs[2:9], px[:, 6], '--', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()
[[0.00268041 0.00154615 0.00139582 0.00280235 0.00280723 0.00142768
  0.00144015]
 [0.00137882 0.00274691 0.00116189 0.00307689 0.00142022 0.0026693
  0.00118179]
 [0.00114741 0.00307596 0.00110675 0.00315649 0.00116595 0.00301502
  0.00112009]
 [0.00109689 0.00315977 0.00108789 0.00317959 0.00110969 0.00310298
  0.00109941]
 [0.00108228 0.0031819  0.00108013 0.00318681 0.00109353 0.00312614
  0.00109107]
 [0.00107715 0.0031881  0.00107662 0.00318933 0.00108795 0.0031326
  0.00108735]
 [0.00107508 0.00319    0.00107495 0.00319031 0.00108575 0.00313456
  0.0010856 ]
 [0.00107417 0.00319065 0.00107414 0.00319073 0.00108478 0.00313522
  0.00108475]
 [0.00107375 0.0031909  0.00107374 0.00319092 0.00108434 0.00313548
  0.00108433]]
Out[8]:
<matplotlib.legend.Legend at 0x1895eda6910>
In [11]:
# Temporal MMS tests
errorst = np.ones((8, 7))
errornorms = np.ones(7)
nx = 100;
nt = 250;
nxs = [25, 50, 100, 200, 400]
nts = [50, 100, 200, 400, 800, 1600, 3200, 6400]
for i in range(8):
    errornorms = numericalcodes(nx, nt)
    errorst[i, 0] = errornorms[0]
    errorst[i, 1] = errornorms[1]
    errorst[i, 2] = errornorms[2]
    errorst[i, 3] = errornorms[3]
    errorst[i, 4] = errornorms[4]
    errorst[i, 5] = errornorms[5]
    errorst[i, 6] = errornorms[6]
    nt = nt*2;
print(errorst)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
# fig.set_size_inches(10.,7.)
plt.plot(nts, np.log(errorst[:, 0]), '-', label = 'CDI')
plt.plot(nts, np.log(errorst[:, 1]), '.-', label = 'CDCN')
plt.plot(nts, np.log(errorst[:, 2]), '-+', label = 'QI')
plt.plot(nts, np.log(errorst[:, 3]), '-*', label = 'QCN')
plt.plot(nts, np.log(errorst[:, 4]), '-x', label = 'eCDI')
plt.plot(nts, np.log(errorst[:, 5]), '--', label = 'eCDCN')
plt.plot(nts, np.log(errorst[:, 6]), '-v', label = 'eQI')
plt.legend(loc = 'best')
plt.xlabel('Grid Numbers')
plt.ylabel('Log Error')

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
pt = np.zeros((6, 7))
for g in range(6):
    for k in range(7):
        pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nts[g]/nts[g+1]))

plt.plot(nts[1:7], pt[:, 0], '-', label = 'CDI')
plt.plot(nts[1:7], pt[:, 1], '.-', label = 'CDCN')
plt.plot(nts[1:7], pt[:, 2], '-+', label = 'QI')
plt.plot(nts[1:7], pt[:, 3], '-*', label = 'QCN')
plt.plot(nts[1:7], pt[:, 4], '-x', label = 'eCDI')
plt.plot(nts[1:7], pt[:, 5], '-v', label = 'eCDCN')
plt.plot(nts[1:7], pt[:, 6], '--', label = 'eQI')
plt.xlabel('Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()
[[4.80797305e-04 1.15463686e-03 4.40796293e-04 1.23423173e-03
  4.92591522e-04 1.12685979e-03 4.47041682e-04]
 [2.70027692e-04 5.23881067e-04 2.23849185e-04 6.02638800e-04
  2.79732016e-04 5.07387787e-04 2.27857313e-04]
 [1.74819852e-04 2.11660470e-04 1.17467654e-04 2.88579416e-04
  1.83282460e-04 2.01049458e-04 1.20428469e-04]
 [1.35616045e-04 6.30726926e-05 6.61680044e-05 1.32099918e-04
  1.43018264e-04 5.72280516e-05 6.86376449e-05]
 [1.20546357e-04 4.42000304e-05 4.26769721e-05 5.43335537e-05
  1.27101641e-04 4.92834621e-05 4.48618322e-05]
 [1.14690533e-04 7.43220864e-05 3.28202700e-05 1.71079475e-05
  1.20701354e-04 8.01449362e-05 3.47624483e-05]
 [1.12262471e-04 9.19383565e-05 2.89506673e-05 1.04211480e-05
  1.17967154e-04 9.75857228e-05 3.06933237e-05]
 [1.11183181e-04 1.01009111e-04 2.74193734e-05 1.73996971e-05
  1.16726232e-04 1.06528305e-04 2.90304857e-05]]
<ipython-input-11-e7ef7324d6f6>:40: RuntimeWarning: invalid value encountered in log
  pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nts[g]/nts[g+1]))
Out[11]:
<matplotlib.legend.Legend at 0x1895ee497c0>
In [12]:
# Printing the spatial grid tests
import pandas as pd
zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Spatial25.csv', index=False)
In [13]:
# Printing the temporal grid tests
zipped2 = list(zip(nts, errorst[:, 0], errorst[:, 1], errorst[:, 2], errorst[:, 3], errorst[:, 4], errorst[:, 5], errorst[:, 6],  pt[:, 0], pt[:, 1], pt[:, 2], pt[:, 3], pt[:, 4], pt[:, 5], pt[:, 6]))

df2 = pd.DataFrame(zipped2, columns=['Temporal Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Temporal250.csv', index=False)
In [15]:
# Dirichlet both sides for spatial 

# Spatial MMS tests
errorsx = np.ones((9, 7))
errornorms = np.ones(7)
nx = 25;
nt = 2000;
nxs = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400]
nts = [250, 500, 1000, 2000, 4000]
for i in range(9):
    errornorms = numericalcodes_dirichlet(nx, nt)
    errorsx[i, 0] = errornorms[0]
    errorsx[i, 1] = errornorms[1]
    errorsx[i, 2] = errornorms[2]
    errorsx[i, 3] = errornorms[3]
    errorsx[i, 4] = errornorms[4]
    errorsx[i, 5] = errornorms[5]
    errorsx[i, 6] = errornorms[6]
    nx = nx*2;
print(errorsx)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 12})
plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI')
plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN')
plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI')
plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN')
plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI')
plt.plot(nxs, np.log(errorsx[:, 5]), '-v', label = 'eCDCN')
plt.plot(nxs, np.log(errorsx[:, 6]), '--', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel(' Log Error (MSE)')
plt.legend()

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 12})
px = np.zeros((7, 7))
for g in range(7):
    for k in range(7):
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))

plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN')
plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI')
plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN')
plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI')
plt.plot(nxs[1:8], px[:, 5], '-v', label = 'eCDCN')
plt.plot(nxs[1:8], px[:, 6], '--', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()

# Printing the spatial grid tests
import pandas as pd
zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Spatial25_Dirichlet.csv', index=False)
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-15-78b68e555be1> in <module>
      9 nts = [250, 500, 1000, 2000, 4000]
     10 for i in range(9):
---> 11     errornorms = numericalcodes_dirichlet(nx, nt)
     12     errorsx[i, 0] = errornorms[0]
     13     errorsx[i, 1] = errornorms[1]

<ipython-input-4-d04c1cfcce3c> in numericalcodes_dirichlet(nx, nt)
    586 
    587 
--> 588         u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
    589         #print(u)
    590         #print(str(n*dt)+" "+str(max(u)))

<ipython-input-4-d04c1cfcce3c> in numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r)
    365 
    366         # print(A)
--> 367         ke = solve(A, bi)
    368         u_eQI = ((1-np.exp(-ke))/frac)
    369         return u_eQI

~\anaconda3\lib\site-packages\scipy\linalg\basic.py in solve(a, b, sym_pos, lower, overwrite_a, overwrite_b, debug, check_finite, assume_a, transposed)
    205         norm = '1'
    206 
--> 207     anorm = lange(norm, a1)
    208 
    209     # Generalized case 'gesv'

MemoryError: Unable to allocate 1.22 MiB for an array with shape (400, 400) and data type float64
In [ ]:
# Spatial MMS tests 2
errorsx = np.ones((9, 7))
errornorms = np.ones(7)
nx = 20;
nt = 1000;
nxs = [20, 40, 80, 160, 320, 640, 1280, 2560, 5120]
nts = [250, 500, 1000, 2000, 4000]
for i in range(9):
    errornorms = numericalcodes(nx, nt)
    errorsx[i, 0] = errornorms[0]
    errorsx[i, 1] = errornorms[1]
    errorsx[i, 2] = errornorms[2]
    errorsx[i, 3] = errornorms[3]
    errorsx[i, 4] = errornorms[4]
    errorsx[i, 5] = errornorms[5]
    errorsx[i, 6] = errornorms[6]
    nx = nx*2;
print(errorsx)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI')
plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN')
plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI')
plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN')
plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI')
plt.plot(nxs, np.log(errorsx[:, 5]), '-v', label = 'eCDCN')
plt.plot(nxs, np.log(errorsx[:, 6]), '--', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Log Error (MSE)')
plt.legend()

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
px = np.zeros((7, 7))
for g in range(7):
    for k in range(7):
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))

plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN')
plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI')
plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN')
plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI')
plt.plot(nxs[1:8], px[:, 5], '-v', label = 'eCDCN')
plt.plot(nxs[1:8], px[:, 6], '--', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()

# Printing the spatial grid tests
zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Spatial20.csv', index=False)

The initial and boundary condition terms: I.C.: Cisin(pix) Boundary Conditions: Periodic BC $$ U(x, 0) = C_o\cos( \pi(2x -1))$$ $$ U(0, t) = -C_o $$ $$ U(1, t) = -C_o $$

$$ \frac {dU}{dx} (0, t) = -2\pi C_o\sin( \pi(2*0 -1)) e^{-t} = 0$$

$$ \frac {dU}{dx} (1, t) = -2\pi C_o\sin(\pi(2*1 -1)) e^{-t} = 0$$

The derivative terms are:

$$ \frac {dU}{dt} = -C_o\cos(\pi (2*x-1)) e^{-t}$$$$ \frac {dU}{dx} = -C_o 2 \pi \sin(\pi (2*x-1)) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o 4 \pi^2 \cos(\pi (2*x-1)) e^{-t}$$

The new source term will be:

$$ d = \frac {dU}{dt} + v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2}$$$$ d = -C_o\cos(\pi (2*x-1)) e^{-t} - v*2*C_o\pi \sin(\pi (2*x-1)) e^{-t} + D*C_o*4*\pi^2 \cos(\pi (2*x-1)) e^{-t}$$$$ d = -C_o e^{-t}(\cos(\pi (2*x-1)) + v*2*\pi \sin(\pi (2*x-1)) - D*4*\pi^2 \cos(\pi (2*x-1))) $$

For the equation with zero and first order decay terms:

Other details remain same, only d will change as shown below:

$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x} + c*U + do $$$$ IC: U(x, 0) = C_i $$$$ BC - Left: U(0, t) = C_o $$$$ BC - Right: \frac {dU(1, t)}{dx} = 0 $$

$$ The \ manufactured \ solution \ for \ MMS \ are: $$

$$ U(x,t) = C_o \cos(\pi x) e^{-t}$$$$ \frac {dU}{dt} = -C_o\cos(\pi x) e^{-t}$$$$ \frac {dU}{dx} = -C_o\pi \sin(\pi x) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o\pi^2 \cos(\pi x) e^{-t}$$$$ The \ initial \ and \ boundary \ conditions \ for \ MMS \ are:$$$$ U(x, 0) = C_o\cos(\pi x)$$$$ U(0, t) = C_oe^{-t}$$$$ \frac {dU}{dx} (1, t) = -C_o\sin(\pi 1) e^{-t} = 0$$

$$ The \ new \ form \ of \ d \ is \ as: $$$$ d = \frac {dU}{dt} + v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2} - c*U - do $$$$ d = -C_o\cos(\pi x) e^{-t} - v*C_o\pi \sin(\pi x) e^{-t} + D*C_o\pi^2 \cos(\pi x) e^{-t} - c*(C_o \cos(\pi x) e^{-t}) - do $$$$ d = -C_o e^{-t}(\cos(\pi x) + v*\pi \sin(\pi x) - D*\pi^2 \cos(\pi x) + c*\cos(\pi x)) - do $$
In [ ]:
# Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. 
def numericalcodes_complete(nx, nt): 
    # Define the grid
    a = 0.1;
    b = 0.001;
    c = 0.1;
    Co = 2.0;
    L = 1.0
    T = 1.0
    do = 0.0
    dm = np.zeros((nt, nx));
    U_exact = np.zeros((nt, nx));
    # Define x & t
    x = np.linspace(0, L, nx)
    t = np.linspace(0, T, nt)
    for i in range(nt):
        for j in range (nx):
            dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - a*np.pi*np.sin(np.pi*x[j]) - b*(np.pi**2)*np.cos(np.pi*x[j]) + c*np.cos(np.pi*x[j])) - do
    
    # Dirichlet Boundary Condition
    left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    right = np.ones(nt)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    for i in range(nt):
        for j in range (nx):
            U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i])
    fig = plt.figure()  
    plt.imshow(np.transpose(U_exact), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
    plt.colorbar()
    plt.title("Exact")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    # Discretization of space and time
    dx = L/(nx-1)
    dt = T/(nt-1)
    u = Co*np.cos(np.pi*x) #np.zeros(nx);
    def numerical_scheme_CDI(u_CDI, a, b, c, d,  dx, dt, l, r):
        alpha = a*dt/dx
        beta = b*dt/dx**2
        A = np.zeros((len(u_CDI), len(u_CDI)))
        bi = d*dt+u_CDI;

        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta-alpha/2) - c*dt/2
            A[j, j] = 1+2*beta
            A[j, j+1] = -(beta+alpha/2) - c*dt/2

        N = len(u_CDI)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0


        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r


        u_CDI = solve(A, bi)

        return u_CDI

    def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r):

        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDI*frac)
        N = len(u_eCDI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((len(u_eCDI), len(u_eCDI)))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b
        #plt.legend()
        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2));
        ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2)

        de = np.zeros(N+1);

        for xt in range(0, (N+1)):
            if(xt==N):
                de[xt] = ke[xt];
            else:
                de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt])

        bi = de;

        for j in range(1, len(u_eCDI)-1):
            A[j, j-1] = ce[j]
            A[j, j] = 1+ae[j]
            A[j, j+1] = be[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0


        bi[0] = le
        # bi[N] = re


        ke = solve(A, bi)
        u_eCDI = ((1-np.exp(-ke))/frac)
        return u_eCDI

    def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r):

        # frac = 0.1


        # ke = -np.log(1-u_eCDCN*frac)
        # N = len(u_eCDCN)-1

        # A = np.zeros((len(u_eCDCN), len(u_eCDCN)))
        # udterm = -c
        # ydterm = d
        # vdterm = -a
        # wdterm = b


        # uk = udterm*dt*np.ones(N+1); # f
        # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        # vk = -vdterm*dt*np.ones(N+1); # a
        # v2k = -wdterm*dt*np.ones(N+1); #d
        # wk = wdterm*dt*np.ones(N+1); # b


        # v2t = np.zeros(N+1)

        # v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2));
        # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2))

        # dcn = np.zeros(N+1);
        # for t in range(0, N+1):
        #     if(t==(N)):
        #         dcn[t] = ke[t];
        #     elif(t==0):
        #         dcn[t] = ke[t];
        #     else:
        #         dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDCN*frac)
        N = len(u_eCDCN)-1

        A = np.zeros((len(u_eCDCN), len(u_eCDCN)))

        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) 

        bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t)
        acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4)
        ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) -  ((b*dt)/(8*(dx**2)))*(v2t)

        dcn = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==(N)):
                dcn[xt] = ke[xt];
            elif(xt==0):
                dcn[xt] = ke[xt];
            else:
                dcn[xt] = ke[xt] - c*dt + ((d[xt]*frac + c)*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2))
                #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        bi = dcn;

        for j in range(1, N):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = ccn[j]
            A[j, j] = 1+acn[j]
            A[j, j+1] = bcn[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = le
        # bi[N] = re

        # print(A)
        ke = solve(A, bi)
        u_eCDCN = ((1-np.exp(-ke))/frac)
        return u_eCDCN

    def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r):
        alpha = a*dt/(4*dx)
        beta = b*dt/(2*dx**2)

        A = np.zeros((len(u_CDCN), len(u_CDCN)))

        bi = np.zeros(len(u_CDCN))
        for j in range(1, len(u_CDCN)-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = -(beta - alpha)
            A[j, j] = 1+(2*beta - c*dt/2)
            A[j, j+1] = -(beta + alpha)

            bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta - c*dt/2))+u_CDCN[j+1]*(alpha+beta)+u_CDCN[j-1]*(beta-alpha)

        N = len(u)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_CDCN = solve(A, bi)

        return u_CDCN

    def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r):
        frac = 0.1
        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eQI*frac)
        N = len(u_eQI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((N+1, N+1))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b

        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        # yk = (bia*frac)*dt;#c+e+uk   *np.ones((1, int(N/2)+1)) #d
        v2t = np.zeros(N+1)
        # print(np.shape(v2t))
        v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2])

        v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3])
        # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)])
        # print(np.shape(v2t))
        #v2t = np.concatenate((v2t, v2t[0]))
        #print(np.shape(v2t))


        cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2)
        bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2)
        aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2)
        fq = vk/(8*dx) + v2k*v2t/((8*dx)**2)

        # USING LAST node Q with other side variation
        cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2)
        bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2)
        anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2)
        fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2)

        dq = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==N):
                dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2))
            else:
                dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2]


        bi = dq;

        for j in range(1, N-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cq[j]
            A[j, j] = 1+aq[j]
            A[j, j+1] = bq[j]
            A[j, j+2] = fq[j]


        A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anq[N-1]
        A[N-1, N] = bnq[N-1]
        A[N-1, N-3] = fnq[N-1]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = le
        # bi[N] = re


        # print(A)
        ke = solve(A, bi)
        u_eQI = ((1-np.exp(-ke))/frac)
        return u_eQI

    def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r):
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QI), len(u_QI)))
        bi = d*dt+u_QI;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QI)-2):

            A[j, j-1] = cQ
            A[j, j] = 1+aQ
            A[j, j+1] = bQ
            A[j, j+2] = eQ

        N = len(u_QI)-1

        A[N-1, N-2] = cnQ  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ
        A[N-1, N] = bnQ
        A[N-1, N-3] = enQ


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_QI = solve(A, bi)

        return u_QI


    def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r):
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QCN), len(u_QCN)))
        bi = d*dt+u_QCN;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QCN)-2):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cQ/2
            A[j, j] = 1+aQ/2
            A[j, j+1] = bQ/2
            A[j, j+2] = eQ/2
            bi[j] = bi[j] - cQ*u_QCN[j-1]/2 - aQ*u_QCN[j]/2 - bQ*u_QCN[j+1]/2 - eQ*u_QCN[j+2]/2

        N = len(u_QCN)-1

        A[N-1, N-2] = cnQ/2  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ/2
        A[N-1, N] = bnQ/2
        A[N-1, N-3] = enQ/2
        bi[N-1] = bi[N-1] - cnQ*u_QCN[N-2]/2 - anQ*u_QCN[N-1]/2 - bnQ*u_QCN[N]/2 - enQ*u_QCN[N-3]/2


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_QCN = solve(A, bi)

        return u_QCN


    # u_init = u_exact(x, 1.0)

    # Apply the numerical scheme

    # U_FD = np.zeros((nx, nt))
    U_CDI = np.zeros((nx, nt))
    U_CDCN = np.zeros((nx, nt))
    U_QI = np.zeros((nx, nt))
    U_QCN = np.zeros((nx, nt))
    U_eCDI = np.zeros((nx, nt))
    U_eQI = np.zeros((nx, nt))

    U_eCDCN = np.zeros((nx, nt))


    # error_2norm_FD_CDI = np.zeros(nt)
    error_2norm_CDI_CDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_CDCN_eCDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_QI_eQI = np.zeros(nt)

    error_2norm_CDI_eCDI = np.zeros(nt)
    error_2norm_CDI_QI = np.zeros(nt)
    error_2norm_QI_QCN = np.zeros(nt)
    error_2norm_CDCN_QCN = np.zeros(nt)
    error_2norm_CDI_QCN = np.zeros(nt)


    error_2norm_eCDCN_eCDI = np.zeros(nt)
    error_2norm_eCDCN_eQI = np.zeros(nt)
    error_2norm_eCDI_eQI = np.zeros(nt)

    error_2norm_Exact_QI = np.zeros(nt)
    error_2norm_Exact_eQI = np.zeros(nt)

    error_2norm_Exact_eCDCN = np.zeros(nt)

    error_2norm_Exact_eCDI = np.zeros(nt)
    error_2norm_Exact_CDCN = np.zeros(nt)
    error_2norm_Exact_QCN = np.zeros(nt)
    error_2norm_Exact_CDI = np.zeros(nt)
    #error_2norm_FD_CDI = np.zeros(nt)


    # u_FD = u;
    u_CDI=u;
    u_CDCN = u;
    u_eQI = u;
    u_eCDI = u;

    u_QI=u;
    u_QCN=u;
    u_eCDCN=u;

    for n in range(nt):
        # print(n)
        # u_FD = numerical_scheme(u_FD, dx, dt)
        # #print(u)
        # #print(str(n*dt)+" "+str(max(u)))
        # U_FD[:, n] = u_FD
        u_CDI = numerical_scheme_CDI(u_CDI,a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDI[:, n] = u_CDI

        #error = u_FD - u_CDI
        #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)




        u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDI[:, n] = u_eCDI

        error = u_CDI - u_eCDI
        error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_CDCN = numerical_scheme_CDCN(u_CDCN, a, b, c,  dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDCN[:, n] = u_CDCN
        error = u_CDI - u_CDCN
        error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDCN[:, n] = u_eCDCN
        error = u_CDCN - u_eCDCN
        error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QI = numerical_scheme_QI(u_QI, a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QI[:, n] = u_QI
        error = u_CDCN - u_QI
        error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QI
        error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, a, b, c,  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eQI[:, n] = u_eQI

        error = u_QI - u_eQI
        error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDI - u_eQI
        error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eQI
        error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eCDI
        error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QCN = numerical_scheme_QCN(u_QCN, a, b, c, dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QCN[:, n] = u_QCN
        error = u_CDCN - u_QCN
        error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QCN
        error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_QI - u_QCN
        error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        # Comparing with the exact solution
        error = U_exact[n, :] - u_QI
        error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQI = error_2norm_Exact_QI[n]
        
        error = U_exact[n, :] - u_QCN
        error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQCN = error_2norm_Exact_QCN[n]
        
        error = U_exact[n, :] - u_eCDCN
        error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDCN = error_2norm_Exact_eCDCN[n]

        error = U_exact[n, :] - u_eCDI
        error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDI = error_2norm_Exact_eCDI[n]
        
        error = U_exact[n, :] - u_CDCN
        error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDCN = error_2norm_Exact_CDCN[n]

        error = U_exact[n, :] - u_CDI
        error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDI = error_2norm_Exact_CDI[n]
        
        error = U_exact[n, :] - u_eQI
        error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreQI = error_2norm_Exact_eQI[n]
    # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150)
    errornorms = np.ones(7)
    errornorms[0] = errorCDI
    errornorms[1] = errorCDCN
    errornorms[2] = errorQI
    errornorms[3] = errorQCN
    errornorms[4] = erroreCDI
    errornorms[5] = erroreCDCN
    errornorms[6] = erroreQI
    
    plt.rcParams['font.size'] = 18
    # fig = plt.figure()  
    # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    # plt.colorbar()
    # plt.title("FD")
    # plt.xlabel('Time')
    # plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_FD.png', dpi = 900)
    # plt.show()

    fig = plt.figure()  
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    fig = plt.figure()  
    plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eCDI.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eQI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eQI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_QI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_QCN.png', dpi = 900)
    plt.show()



    fig = plt.figure() 
    #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN')
    #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI')
    #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI')
    #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI')
    #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN')

    #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN')
    #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN')
    #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI')
    #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI')
    # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI')
    # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI')

    plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI')


    plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI')
    plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN')
    plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN')
    
    plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI')
    plt.title("Errors")
    plt.legend()
    plt.xlabel('Iterations in Time')
    plt.ylabel('Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Error_Norm.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI')
    #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN')
    #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI')
    #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI')
    #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN')

    #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI')
    #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN')
    #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN')
    #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI')

    #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI')
    plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI')
    plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI')
    plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN')
    plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN')
    plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI')

    plt.legend()
    plt.title("Log errors")
    plt.xlabel('Iterations in Time')
    plt.ylabel('log-Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    plt.show()


    #fig = plt.figure() 
    #plt.plot(U_exact[:, -1], '-+',label = 'U_exact')
    #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI')
    #plt.legend()
    #plt.title("Comparison of U_exact and U_CDI")
    #plt.xlabel('Length')
    #plt.ylabel('U')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    #plt.show()
    plt.figure()
    # print(np.shape(U_exact))
    # print(np.shape(U_eQI))
    plt.plot(U_eQI[:, nt-1], label='eQI'); plt.plot(U_exact[nt-1, :], label='Exact'); plt.legend(); plt.title("At 0.1 seconds - Pe = 100")
    return errornorms
In [ ]:
# Equation with decay terms 

# Spatial MMS tests
errorsx = np.ones((9, 7))
errornorms = np.ones(7)
nx = 25;
nt = 1000;
nxs = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400]
nts = [250, 500, 1000, 2000, 4000]
for i in range(9):
    errornorms = numericalcodes_complete(nx+1, nt+1)
    errorsx[i, 0] = errornorms[0]
    errorsx[i, 1] = errornorms[1]
    errorsx[i, 2] = errornorms[2]
    errorsx[i, 3] = errornorms[3]
    errorsx[i, 4] = errornorms[4]
    errorsx[i, 5] = errornorms[5]
    errorsx[i, 6] = errornorms[6]
    nx = nx*2;
print(errors)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
plt.plot(nxs, errorsx[:, 0], '-', label = 'CDI')
plt.plot(nxs, errorsx[:, 1], '.-', label = 'CDCN')
plt.plot(nxs, errorsx[:, 2], '-+', label = 'QI')
plt.plot(nxs, errorsx[:, 3], '-*', label = 'QCN')
plt.plot(nxs, errorsx[:, 4], '-x', label = 'eCDI')
plt.plot(nxs, errorsx[:, 5], '-v', label = 'eCDCN')
plt.plot(nxs, errorsx[:, 6], '--', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Log Error (MSE)')
plt.legend()

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
px = np.zeros((7, 7))
for g in range(7):
    for k in range(7):
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))

plt.plot(nxs[1:6], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[1:6], px[:, 1], '-.', label = 'CDCN')
plt.plot(nxs[1:6], px[:, 2], '-+', label = 'QI')
plt.plot(nxs[1:6], px[:, 3], '-*', label = 'QCN')
plt.plot(nxs[1:6], px[:, 4], '-x', label = 'eCDI')
plt.plot(nxs[1:6], px[:, 5], '-v', label = 'eCDCN')
plt.plot(nxs[1:6], px[:, 6], '--', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()

# Printing the spatial grid tests

zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Spatial25_Complete.csv', index=False)
In [ ]:
 
In [ ]:
 
In [ ]:
# Temporal MMS tests
errorst = np.ones((8, 7))
errornorms = np.ones(7)
nx = 100;
nt = 50;
nxs = [25, 50, 100, 200, 400]
nts = [50, 100, 200, 400, 800, 1600, 3200, 6400]
for i in range(5):
    errornorms = numericalcodes_complete(nx, nt)
    errorst[i, 0] = errornorms[0]
    errorst[i, 1] = errornorms[1]
    errorst[i, 2] = errornorms[2]
    errorst[i, 3] = errornorms[3]
    errorst[i, 4] = errornorms[4]
    errorst[i, 5] = errornorms[5]
    errorst[i, 6] = errornorms[6]
    nt = nt*2;
print(errors)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
# fig.set_size_inches(10.,7.)
plt.plot(nts, np.log(errorst[:, 0]), '-', label = 'CDI')
plt.plot(nts, np.log(errorst[:, 1]), '.-', label = 'CDCN')
plt.plot(nts, np.log(errorst[:, 2]), '-+', label = 'QI')
plt.plot(nts, np.log(errorst[:, 3]), '-*', label = 'QCN')
plt.plot(nts, np.log(errorst[:, 4]), '-x', label = 'eCDI')
plt.plot(nts, np.log(errorst[:, 5]), '--', label = 'eCDCN')
plt.plot(nts, np.log(errorst[:, 6]), '-v', label = 'eQI')
plt.legend(loc = 'best')
plt.xlabel('Number of Time-steps')
plt.ylabel('Log Error (MSE)')

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
pt = np.zeros((6, 7))
for g in range(6):
    for k in range(7):
        pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nts[g]/nts[g+1]))

plt.plot(nts[1:5], pt[:, 0], '-', label = 'CDI')
plt.plot(nts[1:5], pt[:, 1], '.-', label = 'CDCN')
plt.plot(nts[1:5], pt[:, 2], '-+', label = 'QI')
plt.plot(nts[1:5], pt[:, 3], '-*', label = 'QCN')
plt.plot(nts[1:5], pt[:, 4], '-x', label = 'eCDI')
plt.plot(nts[1:5], pt[:, 5], '-v', label = 'eCDCN')
plt.plot(nts[1:5], pt[:, 6], '--', label = 'eQI')
plt.xlabel('Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()
# Printing the temporal grid tests
zipped2 = list(zip(nts, errorst[:, 0], errorst[:, 1], errorst[:, 2], errorst[:, 3], errorst[:, 4], errorst[:, 5], errorst[:, 6],  pt[:, 0], pt[:, 1], pt[:, 2], pt[:, 3], pt[:, 4], pt[:, 5], pt[:, 6]))

df2 = pd.DataFrame(zipped2, columns=['Temporal Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Temporal250_complete.csv', index=False)
In [ ]:
# For ordered grid increases
# Equation with decay terms 

# Spatial MMS tests
errorsx = np.ones((9, 7))
errornorms = np.ones(7)
nx = 20;
nt = 1000;
nxs = [20, 40, 80, 160, 320, 640, 1280, 2560, 5120]

# nxs = [25, 50, 100, 200, 400]
nts = [250, 500, 1000, 2000, 4000]
for i in range(10):
    nx = nxs[i];
    errornorms = numericalcodes_complete(int(nx)+1, nt+1)
    errorsx[i, 0] = errornorms[0]
    errorsx[i, 1] = errornorms[1]
    errorsx[i, 2] = errornorms[2]
    errorsx[i, 3] = errornorms[3]
    errorsx[i, 4] = errornorms[4]
    errorsx[i, 5] = errornorms[5]
    errorsx[i, 6] = errornorms[6]
    
print(errors)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 12})
# fig.set_size_inches(10.,7.)
plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI')
plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN')
plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI')
plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN')
plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI')
plt.plot(nxs, np.log(errorsx[:, 5]), '--', label = 'eCDCN')
plt.plot(nxs, np.log(errorsx[:, 6]), '-v', label = 'eQI')
plt.legend(loc = 'best')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Log Error (MSE)')

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 12})
px = np.zeros((7, 7))
for g in range(7):
    for k in range(7):
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))

plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN')
plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI')
plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN')
plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI')
plt.plot(nxs[1:8], px[:, 5], '--', label = 'eCDCN')
plt.plot(nxs[1:8], px[:, 6], '-v', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend(loc = 'best')

# Printing the spatial grid tests

zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Spatial25_Complete_Ordered.csv', index=False)
In [ ]:
px = np.zeros((8, 7))
for g in range(8):
    for k in range(7):
        # print((np.log(nxs[g]/nxs[g+1])))
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))


plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 12})
plt.plot(nxs[0:8], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[0:8], px[:, 1], '-.', label = 'CDCN')
plt.plot(nxs[0:8], px[:, 2], '-+', label = 'QI')
plt.plot(nxs[0:8], px[:, 3], '-*', label = 'QCN')
plt.plot(nxs[0:8], px[:, 4], '-x', label = 'eCDI')
plt.plot(nxs[0:8], px[:, 6], '-v', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()
In [ ]:
np.log(0.5)

Variable Coefficients a, b, c

$$ a = \sin(\pi x)$$$$ b = \cos(\pi x) $$$$ c = (a+b)/2$$
In [406]:
# Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. 
def numericalcodes_variable(nx, nt): 
    # Define the grid

    Co = 2.0;
    L = 1.0
    T = 1.0
    do = 0.0
    dm = np.zeros((nt, nx));
    U_exact = np.zeros((nt, nx));
    # Define x & t
    x = np.linspace(0, L, nx)
    t = np.linspace(0, T, nt)
    am = np.zeros((nt, nx))
    bm = np.zeros((nt, nx))
    cm = np.zeros((nt, nx))
    
    
    
    for i in range(nt):
        for j in range (nx):
            am[i, j] =  0.1+0.1*np.sin(np.pi*x[j])# *np.exp(-t[i]);
            bm[i, j] =  0.001+0.001*np.cos(np.pi*x[j])# *np.exp(-t[i]);
            cm[i, j] = (am[i, j]+bm[i, j])/2;
            
            dm[i, j] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*x[j]) - am[i, j]*np.pi*np.sin(np.pi*x[j]) - bm[i, j]*(np.pi**2)*np.cos(np.pi*x[j]) + cm[i, j]*np.cos(np.pi*x[j])) - do
    
    # Dirichlet Boundary Condition
    left = Co * np.exp(-t); # (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    right = np.ones(nt)# (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))*np.ones(nt)
    for i in range(nt):
        for j in range (nx):
            U_exact[i, j] = Co*np.cos(np.pi*x[j])*np.exp(-t[i])
    fig = plt.figure()  
    plt.imshow(np.transpose(U_exact), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
    plt.colorbar()
    plt.title("Exact")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    # Discretization of space and time
    dx = L/(nx-1)
    dt = T/(nt-1)
    u = Co*np.cos(np.pi*x) #np.zeros(nx);
    def numerical_scheme_CDI(u_CDI, a, b, c, d,  dx, dt, l, r):
        # print(a,b,c,d,l)
        alpha = a*dt/dx
        beta = b*dt/dx**2
        A = np.zeros((len(u_CDI), len(u_CDI)))
        bi = d*dt+u_CDI;

        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta[j]-alpha[j]/2) - c[j]*dt/2
            A[j, j] = 1+2*beta[j]
            A[j, j+1] = -(beta[j]+alpha[j]/2) - c[j]*dt/2

        N = len(u_CDI)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0


        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r


        u_CDI = solve(A, bi)

        return u_CDI

    def numerical_scheme_e_CDI(u_eCDI, dx, dt,a, b, c, d, l, r):

        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDI*frac)
        N = len(u_eCDI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((len(u_eCDI), len(u_eCDI)))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b
        #plt.legend()
        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        ce = -wk/(dx**2) + v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        be = -vk/(2*dx) - wk/(dx**2) - v2t*v2k/(4*(dx**2));
        ae = -u2k*np.exp(ke)/2 + 2*wk/(dx**2)

        de = np.zeros(N+1);

        for xt in range(0, (N+1)):
            if(xt==N):
                de[xt] = ke[xt];
            else:
                de[xt] = ke[xt]*(-u2k[xt]*np.exp(ke[xt])/2+1)+uk[xt]+u2k[xt]*np.exp(ke[xt])

        bi = de;

        for j in range(1, len(u_eCDI)-1):
            A[j, j-1] = ce[j]
            A[j, j] = 1+ae[j]
            A[j, j+1] = be[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0


        bi[0] = le
        # bi[N] = re


        ke = solve(A, bi)
        u_eCDI = ((1-np.exp(-ke))/frac)
        return u_eCDI

    def numerical_scheme_eCDCN(u_eCDCN, dx, dt,a, b, c, d, l, r):

        # frac = 0.1


        # ke = -np.log(1-u_eCDCN*frac)
        # N = len(u_eCDCN)-1

        # A = np.zeros((len(u_eCDCN), len(u_eCDCN)))
        # udterm = -c
        # ydterm = d
        # vdterm = -a
        # wdterm = b


        # uk = udterm*dt*np.ones(N+1); # f
        # u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        # vk = -vdterm*dt*np.ones(N+1); # a
        # v2k = -wdterm*dt*np.ones(N+1); #d
        # wk = wdterm*dt*np.ones(N+1); # b


        # v2t = np.zeros(N+1)

        # v2t[1:N] = (ke[2:N+1]-ke[0:N-1])

        # ccn = -wk/(dx**2)+ v2k*v2t/(4*(dx**2)) + vk/(2*dx) 
        # bcn = -vk/(2*dx)-wk/(dx**2)-v2t*v2k/(4*(dx**2));
        # acn = -u2k*np.exp(ke)/2 + 2*wk/((dx**2))

        # dcn = np.zeros(N+1);
        # for t in range(0, N+1):
        #     if(t==(N)):
        #         dcn[t] = ke[t];
        #     elif(t==0):
        #         dcn[t] = ke[t];
        #     else:
        #         dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        frac = 0.1

        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eCDCN*frac)
        N = len(u_eCDCN)-1

        A = np.zeros((len(u_eCDCN), len(u_eCDCN)))

        v2t = np.zeros(N+1)

        v2t[1:N] = (ke[2:N+1]-ke[0:N-1]) 

        bcn = ((-a*dt)/(4*dx)) - ((b*dt)/(2*(dx**2))) + ((b*dt)/(8*(dx**2)))*(v2t)
        acn = ((b*dt)/(dx**2)) - ((d*frac + c)*dt*np.exp(ke)/4)
        ccn = (a*dt)/(4*dx) - ((b*dt)/(2*(dx**2))) -  ((b*dt)/(8*(dx**2)))*(v2t)

        dcn = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==(N)):
                dcn[xt] = ke[xt];
            elif(xt==0):
                dcn[xt] = ke[xt];
            else:
                dcn[xt] = ke[xt] - c[xt]*dt + ((d[xt]*frac + c[xt])*dt*np.exp(ke[xt]))*(1-ke[xt]/4)+ (1/2)*(((a[xt]*dt)/(2*dx))*(ke[xt+1] - ke[xt-1]) + ((b[xt]*dt)/(dx**2))*(ke[xt+1]+ke[xt-1]-2*ke[xt]) - ((b[xt]*dt)/(4*(dx**2)))*((ke[xt+1] - ke[xt-1])**2))
                #dcn[t] = (ke[t])*((-u2k[t]*np.exp(ke[t]))/4+1)+(uk[t])/2+(u2k[t]*np.exp(ke[t]))/2 + (1/2)*(u2k[t]*np.exp(ke[t]) + uk[t])+(ke[t+1]/2)*(vk[t]/(2*dx)+wk[t]/(dx**2))+(v2k[t]/2)*(((ke[t+1]-ke[t-1])/(2*dx))**2)+(ke[t]/2)*(((u2k[t]/(2*dx))*np.exp(ke[t])) - (2*wk[t]/(dx**2)))+(ke[t-1]/2)*(-vk[t]/(2*dx) + wk[t]/(dx**2))


        bi = dcn;

        for j in range(1, N):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = ccn[j]
            A[j, j] = 1+acn[j]
            A[j, j+1] = bcn[j]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = le
        # bi[N] = re

        # print(A)
        ke = solve(A, bi)
        u_eCDCN = ((1-np.exp(-ke))/frac)
        return u_eCDCN

    def numerical_scheme_CDCN(u_CDCN, a, b, c, d, dx, dt, l, r):
        alpha = a*dt/(4*dx)
        beta = b*dt/(2*dx**2)

        A = np.zeros((len(u_CDCN), len(u_CDCN)))

        bi = np.zeros(len(u_CDCN))
        for j in range(1, len(u_CDCN)-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = -(beta[j] - alpha[j])
            A[j, j] = 1+(2*beta[j] - c[j]*dt/2)
            A[j, j+1] = -(beta[j] + alpha[j])

            bi[j] = d[j]*dt + u_CDCN[j]*(1-(2*beta[j] - c[j]*dt/2))+u_CDCN[j+1]*(alpha[j]+beta[j])+u_CDCN[j-1]*(beta[j]-alpha[j])

        N = len(u)-1
        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_CDCN = solve(A, bi)

        return u_CDCN

    def numerical_scheme_eQI(u_eQI, dx, dt, a, b, c, d, l, r):
        frac = 0.1
        le = -np.log(1-l*frac)
        re = -np.log(1-r*frac)
        ke = -np.log(1-u_eQI*frac)
        N = len(u_eQI)-1
        #plt.plot(u_eCDI, label='u_eCDI')
        #plt.plot(ke, label='ke')
        A = np.zeros((N+1, N+1))
        udterm = -c
        ydterm = d
        vdterm = -a
        wdterm = b

        uk = udterm*dt*np.ones(N+1); # f
        u2k = (-udterm + ydterm*frac)*dt*np.ones(N+1) # c
        vk = -vdterm*dt*np.ones(N+1); # a
        v2k = -wdterm*dt*np.ones(N+1); #d
        wk = wdterm*dt*np.ones(N+1); # b


        # yk = (bia*frac)*dt;#c+e+uk   *np.ones((1, int(N/2)+1)) #d
        v2t = np.zeros(N+1)
        # print(np.shape(v2t))
        v2t[1:N-1] = (7*ke[2:N]-3*ke[1:N-1] - ke[3:N+1] - 3*ke[0:N-2])

        v2t[N-1] = (3*ke[N]+3*ke[N-1]-7*ke[N-2]+ke[N-3])
        # v2t = (ke[1:int(N/2)+1]-ke[0:int(N/2)])
        # print(np.shape(v2t))
        #v2t = np.concatenate((v2t, v2t[0]))
        #print(np.shape(v2t))


        cq = 3*vk/(8*dx) - wk/(dx**2) + 3*v2k*v2t/((8*dx)**2)
        bq = -7*vk/(8*dx) - wk/(dx**2) - 7*v2k*v2t/((8*dx)**2)
        aq = 3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 + 3*v2k*v2t/((8*dx)**2)
        fq = vk/(8*dx) + v2k*v2t/((8*dx)**2)

        # USING LAST node Q with other side variation
        cnq = 7*vk/(8*dx) - wk/(dx**2) + 7*v2k*v2t/((8*dx)**2)
        bnq = -3*vk/(8*dx) - wk/(dx**2) - 3*v2k*v2t/((8*dx)**2)
        anq = -3*vk/(8*dx) + 2*wk/(dx**2) - u2k*np.exp(ke)/2 - 3*v2k*v2t/((8*dx)**2)
        fnq = -vk/(8*dx) - v2k*v2t/((8*dx)**2)

        dq = np.zeros(N+1);
        for xt in range(0, N+1):
            if(xt==N):
                dq[xt] = (ke[xt])#*((-u2k[0, t]*np.log(ke[t]))/2+1)+(uk[0, t])+(u2k[0, t]*np.exp(ke[t])) ; # currenty CD - I - CN terms - + (1/2)*(u2k[0, t]*np.exp(ke[t]) + 2*uk[0, t])+(ke[t+1]/2)*(vk[0, t]/(2*dx[t])+wk[0, t]/(dx[t]**2))+(v2k[0, t]/2)*(((ke[t+1]-ke[t-1])/2*dx[t])**2)+(ke[t]/2)*(1+ ((u2k[0, t]/(2*dx[t]))*np.exp(ke[t])) - (2*wk[0, t]/(dx[t]**2)))+(ke[t-1]/2)*(-vk[0, t]/(2*dx[t]) + wk[0, t]/(dx[t]**2))
            else:
                dq[xt] = ke[xt]*((-u2k[xt]*np.exp(ke[xt])/2)+1)+uk[xt]+u2k[xt]*np.exp(ke[xt]) #(-c[t]*tau[0, t-1])+(1-b[t])*tau[0, t]-a[t]*tau[0, t+1]-e[t]*tau[0, t+2]


        bi = dq;

        for j in range(1, N-1):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cq[j]
            A[j, j] = 1+aq[j]
            A[j, j+1] = bq[j]
            A[j, j+2] = fq[j]


        A[N-1, N-2] = cnq[N-1] # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anq[N-1]
        A[N-1, N] = bnq[N-1]
        A[N-1, N-3] = fnq[N-1]

        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = le
        # bi[N] = re


        # print(A)
        ke = solve(A, bi)
        u_eQI = ((1-np.exp(-ke))/frac)
        return u_eQI

    def numerical_scheme_QI(u_QI, a, b, c, d, dx, dt, l, r):
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QI), len(u_QI)))
        bi = d*dt+u_QI;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QI)-2):

            A[j, j-1] = cQ[j]
            A[j, j] = 1+aQ[j]
            A[j, j+1] = bQ[j]
            A[j, j+2] = eQ[j]

        N = len(u_QI)-1

        A[N-1, N-2] = cnQ[N-1]  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ[N-1]
        A[N-1, N] = bnQ[N-1]
        A[N-1, N-3] = enQ[N-1]


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_QI = solve(A, bi)

        return u_QI


    def numerical_scheme_QCN(u_QCN, a, b, c, d, dx, dt, l, r):
        
        alpha = a*dt
        beta = b*dt

        A = np.zeros((len(u_QCN), len(u_QCN)))
        bi = d*dt+u_QCN;
        cQ = 3*alpha/(8*dx) - beta/(dx**2)
        bQ = -7*alpha/(8*dx) - beta/(dx**2)
        aQ = 3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        eQ = 1*alpha/(8*dx)

        cnQ = 7*alpha/(8*dx) - beta/(dx**2)
        bnQ = -3*alpha/(8*dx) - beta/(dx**2)
        anQ = -3*alpha/(8*dx) + 2*beta/(dx**2) - c*dt
        enQ = -1*alpha/(8*dx)


        for j in range(1, len(u_QCN)-2):
            # A[j, j-2] = -d[j]*dt
            A[j, j-1] = cQ[j]/2
            A[j, j] = 1+aQ[j]/2
            A[j, j+1] = bQ[j]/2
            A[j, j+2] = eQ[j]/2
            bi[j] = bi[j] - cQ[j]*u_QCN[j-1]/2 - aQ[j]*u_QCN[j]/2 - bQ[j]*u_QCN[j+1]/2 - eQ[j]*u_QCN[j+2]/2

        N = len(u_QCN)-1

        A[N-1, N-2] = cnQ[N-1]/2  # applied the a>0 terms for this one
        A[N-1, N-1] = 1+anQ[N-1]/2
        A[N-1, N] = bnQ[N-1]/2
        A[N-1, N-3] = enQ[N-1]/2
        bi[N-1] = bi[N-1] - cnQ[N-1]*u_QCN[N-2]/2 - anQ[N-1]*u_QCN[N-1]/2 - bnQ[N-1]*u_QCN[N]/2 - enQ[N-1]*u_QCN[N-3]/2


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        A[N, N] = 3
        A[N, N-1] = -4
        A[N, N-2] = 1
        bi[N] = 0

        A[0,0] = 1.0
        # A[N,N] = 1.0

        bi[0] = l
        # bi[N] = r



        u_QCN = solve(A, bi)

        return u_QCN


    # u_init = u_exact(x, 1.0)

    # Apply the numerical scheme

    # U_FD = np.zeros((nx, nt))
    U_CDI = np.zeros((nx, nt))
    U_CDCN = np.zeros((nx, nt))
    U_QI = np.zeros((nx, nt))
    U_QCN = np.zeros((nx, nt))
    U_eCDI = np.zeros((nx, nt))
    U_eQI = np.zeros((nx, nt))

    U_eCDCN = np.zeros((nx, nt))


    # error_2norm_FD_CDI = np.zeros(nt)
    error_2norm_CDI_CDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_CDCN_eCDCN = np.zeros(nt)
    error_2norm_CDCN_QI = np.zeros(nt)
    error_2norm_QI_eQI = np.zeros(nt)

    error_2norm_CDI_eCDI = np.zeros(nt)
    error_2norm_CDI_QI = np.zeros(nt)
    error_2norm_QI_QCN = np.zeros(nt)
    error_2norm_CDCN_QCN = np.zeros(nt)
    error_2norm_CDI_QCN = np.zeros(nt)


    error_2norm_eCDCN_eCDI = np.zeros(nt)
    error_2norm_eCDCN_eQI = np.zeros(nt)
    error_2norm_eCDI_eQI = np.zeros(nt)

    error_2norm_Exact_QI = np.zeros(nt)
    error_2norm_Exact_eQI = np.zeros(nt)

    error_2norm_Exact_eCDCN = np.zeros(nt)

    error_2norm_Exact_eCDI = np.zeros(nt)
    error_2norm_Exact_CDCN = np.zeros(nt)
    error_2norm_Exact_QCN = np.zeros(nt)
    error_2norm_Exact_CDI = np.zeros(nt)
    #error_2norm_FD_CDI = np.zeros(nt)


    # u_FD = u;
    u_CDI=u;
    u_CDCN = u;
    u_eQI = u;
    u_eCDI = u;

    u_QI=u;
    u_QCN=u;
    u_eCDCN=u;

    for n in range(nt):
        # print(n)
        # u_FD = numerical_scheme(u_FD, dx, dt)
        # #print(u)
        # #print(str(n*dt)+" "+str(max(u)))
        # U_FD[:, n] = u_FD
        u_CDI = numerical_scheme_CDI(u_CDI,am[n, :], bm[n, :], cm[n, :], dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDI[:, n] = u_CDI

        #error = u_FD - u_CDI
        #error_2norm_FD_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)




        u_eCDI = numerical_scheme_e_CDI(u_eCDI, dx, dt, am[n, :], bm[n, :], cm[n, :],  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDI[:, n] = u_eCDI

        error = u_CDI - u_eCDI
        error_2norm_CDI_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_CDCN = numerical_scheme_CDCN(u_CDCN, am[n, :], bm[n, :], cm[n, :],  dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_CDCN[:, n] = u_CDCN
        error = u_CDI - u_CDCN
        error_2norm_CDI_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        u_eCDCN = numerical_scheme_eCDCN(u_eCDCN, dx, dt, am[n, :], bm[n, :], cm[n, :],  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eCDCN[:, n] = u_eCDCN
        error = u_CDCN - u_eCDCN
        error_2norm_CDCN_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QI = numerical_scheme_QI(u_QI, am[n, :], bm[n, :], cm[n, :], dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QI[:, n] = u_QI
        error = u_CDCN - u_QI
        error_2norm_CDCN_QI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QI
        error_2norm_CDI_QI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_eQI = numerical_scheme_eQI(u_eQI, dx, dt, am[n, :], bm[n, :], cm[n, :],  dm[n, :], left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_eQI[:, n] = u_eQI

        error = u_QI - u_eQI
        error_2norm_QI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDI - u_eQI
        error_2norm_eCDI_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eQI
        error_2norm_eCDCN_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_eCDCN - u_eCDI
        error_2norm_eCDCN_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)


        u_QCN = numerical_scheme_QCN(u_QCN, am[n, :], bm[n, :], cm[n, :], dm[n, :], dx, dt, left[n], right[n])
        #print(u)
        #print(str(n*dt)+" "+str(max(u)))
        U_QCN[:, n] = u_QCN
        error = u_CDCN - u_QCN
        error_2norm_CDCN_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_CDI - u_QCN
        error_2norm_CDI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        error = u_QI - u_QCN
        error_2norm_QI_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)

        # Comparing with the exact solution
        error = U_exact[n, :] - u_QI
        error_2norm_Exact_QI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQI = error_2norm_Exact_QI[n]
        
        error = U_exact[n, :] - u_QCN
        error_2norm_Exact_QCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorQCN = error_2norm_Exact_QCN[n]
        
        error = U_exact[n, :] - u_eCDCN
        error_2norm_Exact_eCDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDCN = error_2norm_Exact_eCDCN[n]

        error = U_exact[n, :] - u_eCDI
        error_2norm_Exact_eCDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreCDI = error_2norm_Exact_eCDI[n]
        
        error = U_exact[n, :] - u_CDCN
        error_2norm_Exact_CDCN[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDCN = error_2norm_Exact_CDCN[n]

        error = U_exact[n, :] - u_CDI
        error_2norm_Exact_CDI[n] = np.linalg.norm(error)/np.sqrt(nx)
        errorCDI = error_2norm_Exact_CDI[n]
        
        error = U_exact[n, :] - u_eQI
        error_2norm_Exact_eQI[n] = np.linalg.norm(error)/np.sqrt(nx)
        erroreQI = error_2norm_Exact_eQI[n]
    # Errors at time 1s for 7 methods(CDI, CDCN, QI, QCN, eCDI, eCDCN, eQI) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150)
    errornorms = np.ones(7)
    errornorms[0] = errorCDI
    errornorms[1] = errorCDCN
    errornorms[2] = errorQI
    errornorms[3] = errorQCN
    errornorms[4] = erroreCDI
    errornorms[5] = erroreCDCN
    errornorms[6] = erroreQI
    
    plt.rcParams['font.size'] = 18
    # fig = plt.figure()  
    # plt.imshow(U_FD, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    # plt.colorbar()
    # plt.title("FD")
    # plt.xlabel('Time')
    # plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_FD.png', dpi = 900)
    # plt.show()

    fig = plt.figure()  
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()


    fig = plt.figure()  
    plt.imshow(U_eCDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eCDI.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_CDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eCDCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eCDCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDCN.png', dpi = 900)
    plt.show()


    fig = plt.figure() 
    plt.imshow(U_eQI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("eQI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_eQI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QI")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_QI.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    plt.imshow(U_QCN, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("QCN")
    plt.xlabel('Time')
    plt.ylabel('Position')

    # fig.set_size_inches(30.,18.)
    # plt.savefig('U_QCN.png', dpi = 900)
    plt.show()



    fig = plt.figure() 
    #plt.plot(error_2norm_CDI_CDCN, '-+', label = 'CDI_CDCN')
    #plt.plot(error_2norm_CDCN_QI,'-*', label = 'CDCN_QI')
    #plt.plot(error_2norm_CDI_QI, '--' ,label = 'CDI_QI')
    #plt.plot(error_2norm_CDI_eCDI, '+' ,label = 'CDI_eCDI')
    #plt.plot(error_2norm_CDCN_eCDCN, 'x' ,label = 'CDCN_eCDCN')

    #plt.plot(error_2norm_QI_QCN, '-', label = 'QI_QCN')
    #plt.plot(error_2norm_CDCN_QCN, '.', label='CDCN_QCN')
    #plt.plot(error_2norm_QI_eQI, 'o', label='QI_eQI')
    #plt.plot(error_2norm_eCDI_eQI, '.-', label='eCDI_eQI')
    # plt.plot(error_2norm_eCDCN_eQI, '--', label='eCDCN_eQI')
    # plt.plot(error_2norm_eCDCN_eCDI, '.', label='eCDCN_eCDI')

    plt.plot(error_2norm_Exact_QI, '*', label = 'Exact_QI')


    plt.plot(error_2norm_Exact_CDI, '+', label = 'Exact_CDI')
    plt.plot(error_2norm_Exact_CDCN, 'x', label = 'Exact_CDDCN')
    plt.plot(error_2norm_Exact_eCDCN, '.-', label = 'Exact_eCDCN')
    
    plt.plot(error_2norm_Exact_eQI, '-', label = 'Exact_eQI')
    plt.title("Errors")
    plt.legend()
    plt.xlabel('Iterations in Time')
    plt.ylabel('Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Error_Norm.png', dpi = 900)
    plt.show()

    fig = plt.figure() 
    # plt.plot(np.log(error_2norm_FD_CDI), '-+',label = 'FD_CDI')
    #plt.plot(np.log(error_2norm_CDI_CDCN), '-*',label = 'CDI_CDCN')
    #plt.plot(np.log(error_2norm_CDCN_QI), '--',label = 'CDCN_QI')
    #plt.plot(np.log(error_2norm_CDI_eCDI), '-+' ,label = 'CDI_eCDI')
    #plt.plot(np.log(error_2norm_CDCN_eCDCN), '-x' ,label = 'CDCN_eCDCN')

    #plt.plot(np.log(error_2norm_CDI_QI), '-',label = 'CDI_QI')
    #plt.plot(np.log(error_2norm_QI_QCN), '-.',label = 'QI_QCN')
    #plt.plot(np.log(error_2norm_CDCN_QCN), '.',label = 'CDCN_QCN')
    #plt.plot(np.log(error_2norm_QI_eQI), 'o', label='QI_eQI')

    #plt.plot(np.log(error_2norm_eCDI_eQI), '.-', label='eCDI_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eQI), '--', label='eCDCN_eQI')
    #plt.plot(np.log(error_2norm_eCDCN_eCDI), '.', label='eCDCN_eCDI')
    plt.plot(np.log(error_2norm_Exact_QI), '*', label = 'Exact_QI')
    plt.plot(np.log(error_2norm_Exact_CDI), '+', label = 'Exact_CDI')
    plt.plot(np.log(error_2norm_Exact_CDCN), 'x', label = 'Exact_CDCN')
    plt.plot(np.log(error_2norm_Exact_eCDCN), '.-', label = 'Exact_eCDCN')
    plt.plot(np.log(error_2norm_Exact_eQI), '-', label = 'Exact_eQI')

    plt.legend()
    plt.title("Log errors")
    plt.xlabel('Iterations in Time')
    plt.ylabel('log-Error-2-Norm')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    plt.show()


    #fig = plt.figure() 
    #plt.plot(U_exact[:, -1], '-+',label = 'U_exact')
    #plt.plot(U_CDI[:, -1], '-*',label = 'U_CDI')
    #plt.legend()
    #plt.title("Comparison of U_exact and U_CDI")
    #plt.xlabel('Length')
    #plt.ylabel('U')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('Log-ErrorNorm.png', dpi = 900)
    #plt.show()
    plt.figure()
    # print(np.shape(U_exact))
    # print(np.shape(U_eQI))
    plt.plot(U_eQI[:, nt-1], label='eQI'); plt.plot(U_exact[nt-1, :], label='Exact'); plt.legend(); plt.title("At 0.1 seconds - Pe = 100")
    return errornorms
In [407]:
# For ordered grid increases with variable coefficients with x and t
# Equation with first order decay term 

# Spatial MMS tests
errorsx = np.ones((9, 7))
errornorms = np.ones(7)
nx = 20;
nt = 1000;
#nxs = np.linspace(20, 400, 10)

nxs = [20, 40, 80, 160, 320, 640, 1280, 2560, 5120]
nts = [250, 500, 1000, 2000, 4000]
for i in range(10):
    nx = nxs[i]
    errornorms = numericalcodes_variable(int(nx)+1, nt+1)
    errorsx[i, 0] = errornorms[0]
    errorsx[i, 1] = errornorms[1]
    errorsx[i, 2] = errornorms[2]
    errorsx[i, 3] = errornorms[3]
    errorsx[i, 4] = errornorms[4]
    errorsx[i, 5] = errornorms[5]
    errorsx[i, 6] = errornorms[6]
    
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 12})
# fig.set_size_inches(10.,7.)
plt.plot(nxs, np.log(errorsx[:, 0]), '-', label = 'CDI')
plt.plot(nxs, np.log(errorsx[:, 1]), '.-', label = 'CDCN')
plt.plot(nxs, np.log(errorsx[:, 2]), '-+', label = 'QI')
plt.plot(nxs, np.log(errorsx[:, 3]), '-*', label = 'QCN')
plt.plot(nxs, np.log(errorsx[:, 4]), '-x', label = 'eCDI')
plt.plot(nxs, np.log(errorsx[:, 5]), '--', label = 'eCDCN')
plt.plot(nxs, np.log(errorsx[:, 6]), '-v', label = 'eQI')
plt.legend(loc = 'best')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Log Error (MSE)')

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 12})
px = np.zeros((7, 7))
for g in range(7):
    for k in range(7):
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))

plt.plot(nxs[1:8], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[1:8], px[:, 1], '-.', label = 'CDCN')
plt.plot(nxs[1:8], px[:, 2], '-+', label = 'QI')
plt.plot(nxs[1:8], px[:, 3], '-*', label = 'QCN')
plt.plot(nxs[1:8], px[:, 4], '-x', label = 'eCDI')
plt.plot(nxs[1:8], px[:, 5], '--', label = 'eCDCN')
plt.plot(nxs[1:8], px[:, 6], '-v', label = 'eQI')
plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend(loc = 'best')

# Printing the spatial grid tests

zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], errorsx[:, 2], errorsx[:, 3], errorsx[:, 4], errorsx[:, 5], errorsx[:, 6], px[:, 0], px[:, 1], px[:, 2], px[:, 3], px[:, 4], px[:, 5], px[:, 6]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Spatial25_Variable_Ordered.csv', index=False)
In [409]:
# Temporal MMS tests
errorst = np.ones((9, 7))
errornorms = np.ones(7)
nx = 100;
nt = 50;
nxs = [25, 50, 100, 200, 400]
nts = [25, 50, 100, 200, 400, 800, 1600, 3200, 6400]
for i in range(7):
    nt = nts[i]
    errornorms = numericalcodes_complete(nx, nt)
    errorst[i, 0] = errornorms[0]
    errorst[i, 1] = errornorms[1]
    errorst[i, 2] = errornorms[2]
    errorst[i, 3] = errornorms[3]
    errorst[i, 4] = errornorms[4]
    errorst[i, 5] = errornorms[5]
    errorst[i, 6] = errornorms[6]
    
print(errors)
plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
# fig.set_size_inches(10.,7.)
plt.plot(nts, np.log(errorst[:, 0]), '-', label = 'CDI')
plt.plot(nts, np.log(errorst[:, 1]), '.-', label = 'CDCN')
plt.plot(nts, np.log(errorst[:, 2]), '-+', label = 'QI')
plt.plot(nts, np.log(errorst[:, 3]), '-*', label = 'QCN')
plt.plot(nts, np.log(errorst[:, 4]), '-x', label = 'eCDI')
plt.plot(nts, np.log(errorst[:, 5]), '--', label = 'eCDCN')
plt.plot(nts, np.log(errorst[:, 6]), '-v', label = 'eQI')
plt.legend(loc = 'best')
plt.xlabel('Grid Numbers')
plt.ylabel('Log Error(MSE)')

plt.figure()
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 15})
# fig.set_size_inches(10.,7.)
pt = np.zeros((7, 7))
for g in range(7):
    for k in range(7):
        pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nxs[g]/nxs[g+1]))

plt.plot(nts[1:8], pt[:, 0], '-', label = 'CDI')
plt.plot(nts[1:8], pt[:, 1], '.-', label = 'CDCN')
plt.plot(nts[1:8], pt[:, 2], '-+', label = 'QI')
plt.plot(nts[1:8], pt[:, 3], '-*', label = 'QCN')
plt.plot(nts[1:8], pt[:, 4], '-x', label = 'eCDI')
plt.plot(nts[1:8], pt[:, 5], '-v', label = 'eCDCN')
plt.plot(nts[1:8], pt[:, 6], '--', label = 'eQI')
plt.xlabel('Number of Time-Steps')
plt.ylabel('Observed Order')
plt.legend(loc = 'best')
# Printing the temporal grid tests
zipped2 = list(zip(nts, errorst[:, 0], errorst[:, 1], errorst[:, 2], errorst[:, 3], errorst[:, 4], errorst[:, 5], errorst[:, 6],  pt[:, 0], pt[:, 1], pt[:, 2], pt[:, 3], pt[:, 4], pt[:, 5], pt[:, 6]))

df2 = pd.DataFrame(zipped2, columns=['Temporal Discretization Numbers', 'CDI_Error', 'CDCN_Error', 'QI_Error', 'QCN_Error', 'eCDI_Error', 'eCDCN_Error', 'eQI_Error', 'CDI_Order', 'CDCN_Order', 'QI_Order', 'QCN_Order', 'eCDI_Order', 'eCDCN_Order', 'eQI_Order'])

df2.to_csv('Analysis_Temporal250_complete.csv', index=False)
[[4.80797305e-04 1.15463686e-03 4.40796293e-04 1.23423173e-03
  4.92591522e-04 1.12685979e-03 4.47041682e-04]
 [2.70027692e-04 5.23881067e-04 2.23849185e-04 6.02638800e-04
  2.79732016e-04 5.07387787e-04 2.27857313e-04]
 [1.74819852e-04 2.11660470e-04 1.17467654e-04 2.88579416e-04
  1.83282460e-04 2.01049458e-04 1.20428469e-04]
 [1.35616045e-04 6.30726926e-05 6.61680044e-05 1.32099918e-04
  1.43018264e-04 5.72280516e-05 6.86376449e-05]
 [1.20546357e-04 4.42000304e-05 4.26769721e-05 5.43335537e-05
  1.27101641e-04 4.92834621e-05 4.48618322e-05]]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-409-9df30d52c836> in <module>
     34 for g in range(5):
     35     for k in range(7):
---> 36         pt[g, k] = np.log((np.abs(errorst[g+2, k]) - np.abs(errorst[g+1, k]))/(np.abs(errorst[g+1, k]) - np.abs(errorst[g, k])))/(np.log(nxs[g]/nxs[g+1]))
     37 
     38 plt.plot(nxs[0:5], pt[:, 0], '-', label = 'CDI')

IndexError: list index out of range
<Figure size 432x288 with 0 Axes>
In [ ]:
 
In [ ]: